Source: Wikimedia Commons.
The weights wi are determined by imposing conditions. First, we express that the estimate Ψ(r) must be unbiased (mean of the difference between the true function and the interpolant equal to zero)
(4.410)
where m is the mean of F(r). This implies that the sum of the weights is equal to 1:
(4.411)
Second, we minimize the mean square error between the true and interpolant functions, subject to the unbiased condition, by minimizing
(4.412)
where μ is a so-called Lagrange multiplier (a parameter used in optimization theory to find the local maxima and minima of a function subject to equality constraints). This is achieved by equating to zero the partial derivative of G with respect to wi and μ. After some algebraic manipulations, one obtains the following linear system
(4.413)
or, in matrix form,
(4.414)
Here, index 0 refers to the unsampled position r where we seek an estimate of the interpolant Ψ. The solution of the system provides the value for an unsampled point r of the optimal weights wi and the Lagrangian multiplier μ. The variance for ordinary kriging
(4.415)
is a measure of the interpolation error.
The kriging method has the advantage of providing an approximation of the true function that is based on a spatial statistical analysis of the data, and therefore automatically accounts for the possible clustering between data points. Its strength stems from the use of the semi-variance γ rather than geometric distances. It is particularly suited for situations in which data are sparse. A drawback of the method is its complexity and computational burden. Finally, the assumption of stationarity is not always valid. The ordinary kriging method is illustrated by Figure 4.33b.
4.16.5 Correction for Local Effects
Some corrections in the interpolation process may be introduced to account for local effects. We examine here two different approaches that address this issue, taking as an example the interpolation of chemical concentrations measured by a network of monitoring stations.
Innovation kriging (IK) method
The mapping of chemical concentrations by the ordinary kriging method described in the previous section accounts exclusively for observations at monitoring stations. The interpolation can potentially be improved, and specifically account for small-scale patterns such as chemical plumes, by including additional information from a chemical transport model. In the IK method (Blond et al., 2003), this is performed by replacing (4.406) with
(4.416)
where Fi denotes the chemical concentration field observed at N stations located at points ri, and M(r) is a first guess (prior) of the interpolated field provided by the model (background value). The correction term applied to this prior estimate of the field
Hi ≡ H(ri) = Fi − M(ri)
(4.417)
is called the increment or innovation. Again, the weight functions wi need to be determined by optimization. If the spatial resolution of the model is sufficiently high, the prior estimate M(r) may exhibit detailed patterns not seen by the monitoring stations.
In the IK technique, the weighting functions wi are determined by applying the ordinary kriging procedure to the innovation Hi rather than to the observations Fi. The system resulting from the optimization process is similar to system (4.413), but with the covariances applying to the innovation Hi rather than to the observed concentrations Fi. Figure 4.34 illustrates how the IK method can outperform the ordinary kriging technique by providing model information on small-scale features undamped by the observations.
Figure 4.34 Production of a map of surface ozone mixing ratios [ppbv] in the vicinity of Paris based on surface observations at several monitoring stations on July 17, 1999. The measured mixing ratios at these stations are indicated next to the black dots. (a): Analysis by the ordinary kriging method based on the information provided only by the monitoring stations. (b): Analysis by the innovation kriging method in which data from the observing stations are combined with prior information from a chemical transport model. The model predicts the presence of an ozone plume downwind of Paris with a maximum concentration of 115 ppbv. The maps produced by the two techniques are very different, which highlights in this case the importance of the prior information provided by the model.
From Blond et al. (2003).
Local empirical corrections
Another approach to account for local influences missing from the original interpolation process is to use additional statistical information related to some identified forcing factor. Consider for example the mapping of surface ozone pollution in an urban area. Ozone may be locally titrated by conversion to NO2 near emission hotspots of nitric oxide (NO) following reaction (3.18). This titration is important to characterize for population exposure but occurs at scales ~1 km that cannot be properly captured by the ozone observation network. Instead, one can establish an empirical statistical relationship between ozone concentrations and some relevant parameter for which finer mapping is available, such as land use or population density (Figure 4.35), and use this relationship to estimate smaller-scale features. The empirical relationship is first used to “detrend” the original observational data (remove the actual contribution of the small-scale forcing factors on the observed concentrations). The detrended values are then interpolated through one of the methods described above. The local influence is then reintroduced by a “retrending” procedure based on the empirical relationship. Results shown in Figure 4.35 for NO2 and ozone highlight the importance of the local empirical correction in densely populated areas.
Figure 4.35 Empirical correction applied to improve air quality estimates in Belgium. (a): Trend functions for NO2 and O3: average maximum one-hour concentration values as a function of land-use parameter β (here for weekday summer values between 2001 and 2006). Factor β represents a weighted and normalized sum of land-use indicators including, for example, urban fabric, industrial areas, road and rail networks, arable land, agricultural areas, forests, wetlands, etc. (b): Annual mean NO2 surface concentration for year 2006 obtained by the ordinary kriging interpolation method (right) and adjusted to account for the local effects of land use (left). The small circles on the map indicate the location of the monitoring stations. (c): Same, but for surface ozone.
From Janssen et al. (2008).
4.16.6 Conservative Remapping
Data available in one coordinate system must often be transformed to a different coordinate system. The transfer process is generally based on the interpolation of the physical quantities from a source grid to a target or destination grid, an operation called remapping or regridding. A simple conservative area-weighted method is to use interpolation weights provided by the fractional overlap between source and destination grid cells (Figure 4.36). We denote by Ai the area of source grid cell i, by Ak the area of target grid cell k, and by Ai,k the overlap between the two. By definition,
(4.418)
where Nk is the number of cells i that fall within the destination grid cell k.
Figure 4.36 Schematic representation of remapping from a source grid (rectangular, in blue) to a destination grid (Mercator projection, in green). The stippled area Ai,k represents the overlap between source grid cell i (area Ai) and destination cell k (area Ak).
If Fi(r) is an intensive variable in a source grid cell i (i.e., a quantity whose value is independent of the size of the grid cell), and r is the spatial coordinate for the source grid, the mean value of the variable in the destination grid cell k is given by
(4.419)
If we assume that the mean value of Fi is the same for the entire source grid cell and for the overlapping area with destination grid cell k, so that
(4.420)
then we can compute simply as
(4.421)
If Fi is an extensive variable (i.e., a quantity whose value is proportional to the grid cell area), then cumulative values must be used in the transform
ation and the mean quantity for the destination grid cell k is given by
(4.422)
The above expressions represent a first-order area-weighted scheme that is generally good enough when a single remapping operation needs to be done. The calculation of overlapping area Ai,k is not always straightforward when the source and destination grids use different geographic projections. This problem can be addressed by dividing the source grid cell into a large number of mini-cells with area Am such that Am ≪ Ak. The mini-cells are attributed to point geographical locations and one just needs to count the number Ni,k of mini-cells that fall within destination cell k to derive
Ai , k ≈ Ni , kAm
(4.423)
The first-order remapping method discussed above can be improved by expanding Fi(r) as a Taylor series around the centroid ri of source cell i
(4.424)
where ∇iF represents the spatial gradient of F in cell i, and
(4.425)
The regridded field is now second-order accurate if ∇iF is at least a first-order approximation of the gradient, as may be derived from the difference between adjacent grid cells. We then obtain
(4.426)
The integral on the right-hand side can be approximated in simple ways. The second-order remapping scheme is useful to reduce error involved in repeated interconversions between a source grid and destination grid. When a large number of back-and-forth remappings between the two grids are performed, the signature of the coarse grid on the fine grid becomes visible in the case of the first-order interpolation. The second-order remapping retains better the shape of the original function.
Remapping the three components of the wind velocities between two atmospheric grids requires special attention. The components of the wind velocities are related to each other by the continuity equation for air, and errors introduced in the interpolation process may therefore lead to the violation of mass conservation. The problem is addressed by multiplying the wind velocity components on the source grid by the local air concentration (or pressure), and by applying the conservative interpolation procedure to the resulting mass flux. The interpolated field on the destination grid is then divided by the air concentration (or pressure) to derive the interpolated wind velocity.
References
Andrews, D. G., Holton, J. R. and Leong, C. B. (1987) Middle Atmosphere Dynamics, Academic Press, Orlando.
Andrews D. G. and McIntyre M. E. (1976) Planetary waves in horizontal and vertical shear: The generalized Eliassen–Palm relation and the mean zonal acceleration, J. Atmos. Sci., 33, 2031–2048.
Andrews D. G. and McIntyre M. E. (1978) An exact theory of nonlinear waves on a Lagrangian-mean flow, J. Fluid Mech., 89, 609–646.
Arakawa A. and Lamb V. R. (1977) Computational design of the basic dynamical processes of the UCLA general circulation model, Methods Comput. Phys., 17, 174–265.
Asselin R. (1972) Frequency filter for time integrations, Mon. Wea. Rev., 100, 487–490.
Banks P. M. and Kockarts G. (1973) Aeronomy, Academic Press, New York.
Bishop C. M. (1995) Neural Networks for Pattern Recognition, Clarendon Press, Oxford.
Blond N., Bel L., and Vautard R. (2003) Three-dimensional ozone data analysis with an air quality model over the Paris area, J. Geophys. Res, 108 (D23), 4744, doi:10.1029/2003JD003679.
Boughton B. A., Delarentis J. M., and Dunn W. W. (1987) A stochastic model of particle dispersion in the atmosphere, Boundary-Layer Meteorology, 80, 147–163.
Bourke W. (1974) A multi-level spectral model: I.Formulation and hemispheric integrations, Mon. Wea. Rev., 102, 687–701.
Boyd J. (1976) The noninteraction of waves with the zonally averaged flow on a spherical earth and the interrelationships of eddy fluxes of energy, heat and momentum, J. Atmos. Sci., 33, 2285–2291.
Boyd J. (1998) Two comments on filtering (artificial viscosity) for Chebyshev and Legendre spectral and spectral element methods: Preserving boundary conditions and interpretation of the filter as a diffusion. J. Comput. Phys., 143, 283–288.
Brasseur G. P. and Solomon S. (2005) Aeronomy of the Middle Atmosphere: Chemistry and Physics of the Stratosphere and Mesosphere, 3rd edition, Springer, Amsterdam.
Canuto V. M., Goldman I., and Hubickyj O. (1984) A formula for the Shakura–Sunyaev turbulent viscosity parameter, Astrophys. J., 280, L55–L88, doi: 10.1086/184269.
Chapman S. and Cowling T. G. (1970) The Mathematical Theory of Non-uniform Gases, 3rd edition, Cambridge University Press, Cambridge.
Charney J. G., Fjörtoft R., and von Neumann J. (1950) Numerical integration of the barotropic vorticity equation, Tellus, 2, 237–254.
Courant R. (1943) Variational methods for the solution of problems of equilibrium and vibrations, Bull. Amer. Math. Soc., 49, 1–23.
Cressie N. (1993) Statistics for Spatial Data, Wiley, Chichester.
Cruse H. (2006) Neural Networks as Cybernetic Systems, 2nd edition, Brains, Mind and Media, Bielefeld.
de Bruyns Kops S. M., Riley J. J., and Kosaly G. (2001) Direct numerical simulation of reacting scalar mixing layers, Phys. Fluids, 13, 5, 1450–1465.
Dietachmayer G. S. and Droegemeier K. K. (1992) Application of continuous dynamic grid adaption techniques to meteorological modelling. Part I: Basic formulation and accuracy, Mon. Wea. Rev., 120, 1675–1706.
Durran D. R. (2010) Numerical Methods for Fluid Dynamics, Springer, Amsterdam.
Eliasen E., Machenhauer B., and Rasmussen E. (1970) On a Numerical Method for Integration of the Hydrodynamical Equations with a Spectral Representation of the Horizontal Fields, Institute of Theoretical Meteorology, University of Copenhagen, Copenhagen.
Enting I. G. (2000) Green’s function methods of tracer inversion, Geophys. Monograph., 114, 19–31,
Fournier A., Taylor M. A., and Tribbia J. (2004) The spectral element atmosphere model (SEAM): High resolution parallel computation and localized resolution of regional dynamics, Mon. Wea. Rev., 132, 726–748.
Fritsch F. N. and Carlson R. E. (1980) Monotone piecewise cubic interpolation, SIAM J. Numer. Anal., 17, 238–246.
Garcia, R. and Solomon S. (1983) A numerical model of the zonally averaged dynamical and chemical structure of the middle atmosphere, J. Geophys. Res. 88, 1379–1400.
Garcia-Menendez F., Yano A., Hu Y., and Odman M. T. (2010) An adaptive grid version of CMAQ for improving the resolution of plumes, Atmos. Pollut. Res., 1, 239–249.
Gentile M., Courbin F., and Meylan G. (2013) Interpolating point spread function anisotropy, Astronomy & Astrophysics, 549, A1.
Hall T. M. and Plumb R. A. (1994) Age of air as a diagnostic of transport, J. Geophys. Res., 99, 1059–1070.
Haltiner G. J. and Williams R. T. (1980) Numerical Prediction and Dynamic Meteorology, Wiley, Chichester.
Hubbard M. E. (2002) Adaptive mesh refinement for three-dimensional off-line tracer advection over the sphere, Int. J. Numer. Methods Fluids, 40, 369–377.
Jablonowski C. and Williamson D. L. (2011) The pros and cons of diffusion, filters and fixers in atmospheric general circulation models. In Numerical Techniques for Global Atmospheric Models (Lauritzen P. H., Jablonowski C., Taylor M. A., and Nair R. D., eds.), Springer-Verlag, Berlin.
Jacob D. J. (1999) Introduction to Atmospheric Chemistry, Princeton University Press, Princeton, NJ.
Janssen S., Dumont G., Fierens F., and Mensink C. (2008) Spatial interpolation of air pollution measurements using CORINE land cover data, Atmos. Env., 42, 4884–4903.
Kasahara A. (1974) Various vertical coordinate systems used for numerical weather prediction, Mon. Wea. Rev., 102, 504–522.
Krige D. G. (1951) A statistical approach to some basic mine valuation problems on the Witwatersrand, J. Chem., Metal. Mining Soc. South Africa, 52, 119–139.
Langevin P. (1908) On the theory of Brownian motion, C. R. Acad. Sci. (Paris), 146, 530–533.
Lanser D. and Verwer J. G. (1998) Analysis of operator splitting for advection–diffusion–reaction problems from air pollution mo
deling. CWI Report MAS-R9805.
Lanser D. and Verwer J. G. (1999) Analysis of operator splitting for advection–diffusion-reaction problems from air pollution modelling, J. Comput. Appl. Math., 111, 210–216.
Laprise R. (1992) The resolution of global spectral models. Bull. Amer. Meteor. Soc., 73, 1453–1454.
Lauritzen P. H., Jablonowski C., Taylor M. A., and Nair R. D. (2011) Numerical Techniques for Global Atmospheric Models, Springer, Amsterdam.
Lin J. C. (2012) Lagrangian modeling of the atmosphere: An introduction. In Lagrangian Modeling of the Atmosphere (Lin J., Brunner D., Gerbig C., et al., eds.), American Meteorological Union, Washington, DC.
Lin J. C., Gerbig C., Wofsy S. C., et al. (2003) A near-field tool for simulating the upstream influence of atmospheric observations: The Stochastic Time-Inverted Lagrangian Transport (STILT) model, J. Geophys. Res., 108(D16), 4493, doi:10.1029/202JD003161.
Lin J. S. and Hildemann L. (1996) Analytical solutions of the atmospheric diffusion equation with multiple sources and height-dependent wind speed and eddy diffusivities, Atmos. Environ., 30(2), 239–254.
Modeling of Atmospheric Chemistry Page 26