Towards Hyper-Dimensional Variography Using the Product-Sum Covariance Model

Modeling hyper-dimensional spatial variability is a complex task from both practical and theoretical standpoints. In this paper we develop a method for modeling hyper-dimensional covariance (variogram) structures using the product-sum covariance model initially developed to model spatio-temporal variability. We show that the product-sum model can be used recursively up to an arbitrarily large number of dimensions while preserving relative modeling simplicity and yielding valid covariance models. The method can be used to model variability in anisotropic conditions with multiple axes of anisotropy or when temporal evolution is involved, and thus is applicable to “full anisotropic 3D+time” situations often encountered in environmental sciences. It requires fewer assumptions than the traditional product-sum modeling approach. The new method also presents an alternative to classical approaches to modeling zonal anisotropy and requires fewer parameters to be estimated from data. We present an example by applying the method in conjunction with ordinary kriging to map photosynthetically-active radiation (PAR) for 2006, in Oklahoma, CA, USA and to explore effects of spatio-temporal variability in PAR on gross primary productivity (GPP).


Introduction
Spatial scaling of ecological processes has been a central theme in ecological modeling.Variograms, defined as the variance of the difference between two random variables of the given random function at two locations over the domain [1,2], represent an important tool having applications to the analysis of spatial variability and spatial scaling.They are commonly modeled by creating a graph relating the variance of the difference in value of a variable at pairs of sample points to the separation distance between those pairs.Variography is usually used in tandem with another geostatistical technique-kriging, to obtain kriged estimates of the field at locations of interest [1].
Several theoretical variogram models are usually utilized for modeling simple, one-dimensional, isotropic variograms, including, for example, exponential, Gaussian, and spherical, for cases where spatial variability can be described using only one variogram and one distance measure [1,3].However, more complex situations are typical in environmental science and ecology, and usually oversimplified by using common lower-dimensional theoretical variogram models.Important examples include anisotropy (where variability depends on direction), and spatio-temporal variability.For example, in Tadić et al. [4] the atmosphere was modeled using omnidirectional isotropic vaiogram models, ignoring vertical/horizontal anisotropy, and in Hammerling at al. [5] and Tadić et al.'s [6] observations spanning multiple days were treated as coincident, ignoring the temporal component in spatio-temporal variability.Extending spatial interpolation techniques to the spatio-temporal domain cannot be accomplished simply by adding another dimension [7][8][9].
The origin of the spatial and temporal variation can differ, and that difference can lead to strong zonal and geometric anisotropic behavior [10].For example, modeling an atmospheric domain within h + ∆h, v + ∆v, and t + ∆t, where h, v, and t represent horizontal, vertical, and temporal coordinates, starting from separate one-dimensional variograms would require three independent one-dimensional variograms, given that the atmosphere often exhibits zonal anisotropic properties [4].Similarly, modeling Earth surface temperature within lat + ∆lat, lon + ∆lon, h + ∆h, and t + ∆t would require four independent one-dimensional variogram models, given that temperature is commonly less variable along longitude than latitude; i.e., zonally anisotropic within the horizontal plane.
In this study we show that a relatively easy tweak of the generalized product-sum model [3,22] can be applied to model hyper-dimensional variograms having an arbitrary number of dimensions.
The key advantages of the proposed approach are two-fold: 1) The resulting models retain a relative computational simplicity characteristic for generalized product-sum model.2) We avoid defining a tolerance required by the originally proposed modeling procedure [22], thus simplifying and making the modeling procedure less arbitrary, to facilitate its application to ecological datasets (see Section 2.2.2).
Apart from modeling spatio-temporal variability, modeling hyperdimensional covariance structures allows modeling spatial variability in cases where multiple axes of anisotropy are present, which is of particular interest to atmospheric sciences.For example, when a plume is emitted from point sources, there are usually three principal axes of anisotropy: along the plume in horizontal direction, perpendicular to the plume in horizontal direction, and vertical direction.In a recent study [23] we used the product-sum model to model zonal anisotropy (with only two principal axes) and provided comparison to the classical models of zonal anisotropy based on nested variogram structures, while emphasizing the comparative simplicity of the product-sum model-based approach.
We envision this study as a contribution towards an applied aspect of modeling covariance (variogram) structures in hyper-space.We provide the code and basic suggestions of how to deploy the proposed modeling tools in two different setups (see Sections 2.2.1 and 2.2.2).
After presenting the model and its modifications in Section 2, we describe in Section 3 how we use this technique in conjunction with a canopy photosynthesis model to estimate gross primary productivity (GPP) using solar irradiation measurements from 118 ground stations in the Oklahoma Mesonet [24,25].In this example, we use three independent one-dimensional variogram models (V (3) ) for horizontal, vertical, and temporal axes, thus representing the simplest hyper-dimensional case above the two-dimensional variogram space originally used to develop the product-sum model.
Two main difficulties arise in modeling spatio-temporal and other hyper-dimensional correlations: (a) how to ensure one has a valid model [21,27], and (b) how to fit the model to data [22].An overview of the available classes of spatio-temporal covariance models is given in Montero et al. [3] and [11].Possible models include metric [28], linear [29], product [30], non-separable [31], and generalized product-sum models [22].Here we adopt the product-sum model, which has recently gained popularity due to several advantages (e.g., sequential modeling, computational feasibility) it offers for environmental applications such as satellite observations [26,[32][33][34].
In what follows, the terms covariance and variogram are used interchangeably due to a simple relationship between them [1]: where C h denotes covariance and γ h variogram.To simplify notation, we will use symbol "V" for variogram and "COV" for covariance, and the number in the superscript next to the symbols will denote their dimensionality: V (1) stands for one-dimensional, V (2) two-dimensional variogram, etc., with corresponding covariance COV (1) , COV (2) , etc. Common spatio-temporal modeling (only horizontal spatial and temporal distance considered) belongs to V (2) class.

Original Product-Sum Model and Modeling Procedure
The product-sum variogram (covariance) model [22,35,36] was initially developed to model spatio-temporal covariance structures.Montero et al. [3] provides an exhaustive definition of the model and provides the reader with several examples.The following class of valid product-sum covariance models was introduced in De Cesare et al. [35], and further developed in De Iaco et al. [22]: where C s and C t are valid spatial and temporal covariance models, respectively.De Iaco et al. [22] proved that for positive definiteness, it is sufficient that a 1 > 0, a 2 ≥ 0 and a 3 ≥ 0. The model in Equation ( 2) corresponds to the spatio-temporal variogram shown in Equation (3).In the original procedure, De Iaco et al. [22] estimated separate V (1) spatial (h s = 0) and temporal (h t = 0) variograms using the data, and then combined these models to obtain the final spatio-temporal variogram model: where γ s,t (h s ,0) and γ s,t (0, h t ) are spatio-temporal variograms for h s = 0 and h t = 0, respectively.Parameter k is estimated from the data, which makes the model easy to apply: where k s C s (0) and k t C t (0) are spatial and temporal sills (variances) obtained in modeling of separate spatial and temporal variograms (both V (1) ).The only condition k has to fulfill to create an admissible covariance model is: The original modeling procedure [22] assumed separate modeling of the spatial and temporal covariance (variograms) and their later unification into a spatio-temporal product-sum model.Apart from the apparent simplicity, such an approach posed one problem.Namely, in cases where sparse data are located on an irregular grid, there are often insufficient data to model separate variograms, thus the authors suggested using a tolerance along each dimension, or, in other words, an arbitrarily predetermined maximal allowed distance between data points in order to consider them collocated or coincident.However, if the tolerance is large, it degrades the quality of the V (1) models, and the existence of tolerance per se makes the procedure somewhat subjective.In what follows, we present a new approach that does not require defining tolerances.

Modeling of the Hyper-Dimensional Variogram Based on the Product-Sum Model
Our extension of the product-sum model is based on three previously unexplored properties of Equation ( 2) that allow for extrapolating the approach to hyper-dimensional cases encountered in environmental science and ecology.First, we notice that the validity of the product-sum covariance model using basic COV (1) class of covariances C s and C t is not dependent on the dimensionality of the basic covariance models, and the same applies to corresponding variograms (Equation ( 3)).The minimum and necessary conditions to assure validity of the product-sum model are: (1) that basic COV (1) building blocks are valid models, and (2) that constants a 1 -a 3 are subjected to constraints mentioned earlier in Equation (2) (De Iaco et al. 2001).Thus, the product-sum model validity holds for any basic covariance(s) dimensionality (equivalent to C s and C t in Equation ( 2)), as long as they represent a valid model (for validity criteria, see Chiles and Delfiner, 2012).Subsequently, a COV (2)  class of covariance modeled in the first step starting from basic COV (1) components can be used as a basic covariance in the subsequent steps, thus increasing the dimensionality of the resulting covariance.
Second, Equation ( 2) is symmetric in terms of C s and C t in the sense that if these terms exchange places in the equation, it will still converge to the same expression given that the constants a 1 -a 3 are data-driven, (note that constants a 1 -a 3 are not directly modeled, but rather implicitly through global and partial sills estimated from the data; please see Equations ( 4)-( 8) in De Iaco et al. [22] for further clarifications).We will later show that fitting the model to data can be done sequentially, mimicking the approach from the original paper by De Iaco et al. [22], or, as we propose, all at once, thus avoiding the need to define tolerances.
Third, Equation ( 4) can be rewritten as: Based on the argument from analogy [37], modeling COV (n) type of covariance starting with COV (n−1) and COV (1) would have k n-1 value: For modeling an n-dimensional variogram (covariance), a set of n−1 values of k 1 -k n−1 would have to be estimated from the data.This expression allows us to extend the model to an arbitrary number of dimensions.

Sequential Hierarchical Modeling
Given the above mentioned properties of the product-sum model, assuming we intend to model COV (3) class of the covariance, we could start by separately modeling one COV (2) class of the covariance using the product-sum model, and a COV (1) class of the covariance using basic one-dimensional model, and then combining them again using the product-sum model into a final COV (3) .The approach is depicted in Scheme 1.
Atmosphere 2018, 9, x FOR PEER REVIEW 4 of 11 minimum and necessary conditions to assure validity of the product-sum model are: (1) that basic COV (1) building blocks are valid models, and (2) that constants a1-a3 are subjected to constraints mentioned earlier in Equation (2) (De Iaco et al. 2001).Thus, the product-sum model validity holds for any basic covariance(s) dimensionality (equivalent to Cs and Ct in Equation ( 2)), as long as they represent a valid model (for validity criteria, see Chiles and Delfiner, 2012).Subsequently, a COV (2)  class of covariance modeled in the first step starting from basic COV (1) components can be used as a basic covariance in the subsequent steps, thus increasing the dimensionality of the resulting covariance.
Second, Equation ( 2) is symmetric in terms of Cs and Ct in the sense that if these terms exchange places in the equation, it will still converge to the same expression given that the constants a1-a3 are data-driven, (note that constants a1-a3 are not directly modeled, but rather implicitly through global and partial sills estimated from the data; please see Equations ( 4)- (8) in De Iaco et al. [22] for further clarifications).We will later show that fitting the model to data can be done sequentially, mimicking the approach from the original paper by De Iaco et al. [22], or, as we propose, all at once, thus avoiding the need to define tolerances.
Third, Equation ( 4) can be rewritten as: Based on the argument from analogy [37], modeling COV (n) type of covariance starting with COV (n−1) and COV (1) would have kn-1 value: For modeling an n-dimensional variogram (covariance), a set of n−1 values of k1-kn−1 would have to be estimated from the data.This expression allows us to extend the model to an arbitrary number of dimensions.

Sequential Hierarchical Modeling
Given the above mentioned properties of the product-sum model, assuming we intend to model COV (3) class of the covariance, we could start by separately modeling one COV (2) class of the covariance using the product-sum model, and a COV (1) class of the covariance using basic onedimensional model, and then combining them again using the product-sum model into a final COV (3) .The approach is depicted in Scheme 1.
This approach, analogous to spatio-temporal modeling in Equation ( 2), starts by selecting the first pair of independent one-dimensional variograms to be combined into a product-sum model.In the case of separately modeling variability in horizontal and vertical directions and time, it requires first creating the product-sum model for the first arbitrarily selected pair of dimensions; for example, horizontal and vertical directions: Second, C h,v (h h ,h v ) (2) and C t (h t ) (1) are combined into a final product-sum model: After substituting Equation ( 8) into Equation ( 9) it yields: + a 2 a 5 C h (h h ) (1) + a 3 a 5 C v (h v ) (1) + a 6 C t (h t ) (1)  (10) Two important conclusions follow from Equation (10).First, any n-dimensional variogram (covariance) can be broken down to one-dimensional components following n−1 modeling steps, each comprising the modeling of one product-sum variogram of the lower dimensionality.Second, Equation (10) shows that the order of modeling lower-dimensional variograms is not important, as it converges to the same model.However, the inferred covariance parameters could be slightly different and dependent on the modeling sequence if tolerances are large, though this can be overcome by modeling at-once, as shown in Section 2.2.2.If we start by modeling first the C h,t (h h ,h t ) in Equation ( 8), the initial values for a 1 -a 6 would be different, yet their products (shown in bold) in Equation ( 10) would be the same, as they ultimately come from data, which again points to the symmetrical properties of the product sum model with respect to its basic covariance components.The resulting covariance is guaranteed to yield a valid model, as the product of constants a 1• a 4 that multiplies the first term in Equation ( 10) is always >0, while all other products of constants are ≥0.The four underlined covariance product terms of the final model represent a Hadamard product [38] of two or more positive definite matrices.According to Schur product theorem, a Hadamard product of two positive definite matrices necessarily gives a positive definite matrix [39].Thus, the resulting model is guaranteed to be valid.
In the original sequential modeling procedure starting from one-dimensional variograms, γ s,t (h s ,0) and γ s,t (0, h t ) were modeled by pre-determining a spatial and temporal tolerance within which a pair of data can still be considered collocated or coincident.However, in practice the required tolerance might need to be very large to yield enough data points, which would lead to inaccuracies in the modeling of both one-dimensional and higher-dimensional variograms, given that lower dimensional variograms are building blocks of higher dimensional variogram (see Section 2.2.2 for further discussion).
Furthermore, correctly representing variability along all axes of anisotropy leads to the reduction in the estimated nugget, if nugget-model variograms are used.For example, in a typical spatial-only variography of satellite retrievals, the variography step is preceded by temporal binning, i.e., data from multiple time periods are aggregated and treated as coincident (e.g., [5,6]).The consequence is that ignored variability along the temporal axis ends up being superposed to the spatial variability as its unexplained portion-nugget.
Similarly, the variability within tolerances imposed by the original sequential product-sum modeling procedure [22] gets superposed to the true variability along the actual modeling axis affecting the unexplained portion of variability.By treating not-exactly-collocated and not-exactly-coincident observations as collocated and coincident, modeler allows some of the spatial variability (within tolerance) to affect the separate temporal variogram, and some of the temporal variability to affect the separate spatial variogram, respectively.In this way, separate spatial and temporal variogram models represent an approximation of the separate variogram models that would be obtained using strictly collocated and strictly coincident observations.This issue was not discussed in the original paper by De Iaco et al. [22], although it affects modeling results for the sequential approach (as the nugget has to be modeled in the first step), but not in modeling the "all at once" approach (Section 2.2.2).Thus, we change the approach, and model all parameters at once.

Modeling "All at Once"
To avoid defining a tolerance, we alter the original procedure by estimating all covariance parameters simultaneously.This simultaneous parameter estimation makes the model more applicable to scattered data and data with variable spatial coverage, as is often the case with airborne observations or satellite data.Defining the tolerance followed naturally from the original sequential modeling procedure proposed in De Iaco et al. [22], starting by modeling separate one-dimensional variograms used in subsequent steps to yield a generalized product-sum model.In the original procedure, the first step assumed computing the sample spatial and temporal variograms corresponding to γ s,t (h s ,0)and γ s,t (0,h t ): where r s and r t were, respectively, the vector lag with spatial tolerance δ s and the lag with temporal tolerance δ t .|N(r s )| and |M(r t )| are the cardinalities of the following sets: The need to define the tolerance follows from the fact that, in order to construct separate spatial variogram based on coincident observations ( γs,t (r s , 0)), or separate temporal variogram based on collocated observations ( γs,t (0, r t )), the observational data set in practice has to be dense, or tolerances have to be relatively large to yield enough data for the construction of variograms.In other words, there are often too few collocated pairs of points that cover a wide range of temporal distances, and coincident points that cover a wide range of spatial distances, which imposes increasing the tolerances.The trade-off is that the variabilities captured by separate variograms become virtually higher, as temporal and spatial variability between pairs of points separated within the tolerance gets superposed on spatial and temporal variability of coincident, and collocated points, respectively.However, if the hyperdimensional variograms are fitted to the raw variograms at once, the need to define tolerances vanish.All variogram parameters are estimated simultaneously, including the ones that correspond to separate spatial and temporal variograms.The fitting procedure exploits all found pairs of points at their exact spatial and temporal distances and fits the optimal parameters at once.
The modeling starts by fitting the V (n) variogram corresponding to Equation (3).into an experimental (or alternatively a raw) variogram, and replacing k using Equation (7).Then the V (n−1) variogram component on the right-hand side of the expression is again substituted using Equation (3) and Equation (7), and maximal dimensionality of the basic variograms is again reduced by one.The procedure is recursively repeated until all variogram terms on the right-hand side have dimensionality one.While the expression itself can be very large, the number of parameters to estimate all at once is relatively small and computationally feasible.For example, to model V (3) variogram using common 3-parameter exponential or Gaussian model, it would be required to simultaneously estimate only 9 parameters (nugget, 3 times V (1) sill and range parameters, and two k-parameters).
For temporally evolving anisotropic full 3D space, only 12 estimated parameters would be required.Generally, the number of parameters for an n-dimensional variogram to be estimated is 3n.
From a practical standpoint, in the sequential modeling approach, modeling constants k 1 -k n−1 play important roles, since the sill (2) -sill (n) values are obtained after the k value is optimized in each step.However, in modeling "all at once," k constants are substituted following Equation (7), and do not appear as important entities over the course of modeling.

Application
As an example of how our method can address spatial scaling problems, we explored sensitivity of gross primary productivity (GPP) to spatial variability in photosynthetically-active radiation (PAR).Clouds influence the quantity and quality of light for GPP [40,41].However, gridded meteorological datasets are often too coarse (~100 km) to resolve mesoscale (~10 km) weather systems (e.g., thunderstorms) that affect PAR.Furthermore, modeled GPP calculated from coarse-resolution PAR can be biased due to nonlinearity in the GPP-PAR relationship.We explored the impact of spatial variability in PAR on GPP, using our kriging method to estimate shortwave radiation from Oklahoma Mesonet stations (PAR was approximated as half the downwelling shortwave radiation; [42].We modeled GPP as a function of PAR (black curve in Figure 1), using a two-leaf canopy model applied to winter wheat in Oklahoma [43].
The spatio-temporal modeling of PAR was done under a moving window setup [44] in such a way that spatial covariance for every 5-min time slice was modeled separately, allowing sampling of the data within +/−30 preceding and following 5-min time periods.Thus, problems resulting from the non-stationarity of the PAR mean due to the diurnal cycle were avoided.We assumed that the V (3)  variogram offers an adequate representation of the variability, and treated time, horizontal and vertical distances as orthogonal components in variogram space, and thus modeled separately.The vertical and temporal directions were modeled using Gaussian, and horizontal using an exponential model.After modeling the covariance, we kriged the PAR to a regular ~10 × 10 km grid mimicking the kriging approach in Tadić et al. [6], and obtained the average PAR for May 2006 shown in Figure 1.
We estimated the impact of spatial variability in PAR on GPP, using GPP calculated from spatially-averaged PAR as a basis for comparison.Due to nutrient or other limitations, GPP can become light-saturated, such that further increases in PAR yield diminishing returns for GPP at high PAR values.This saturation effect is characterized by the nonlinear shape of the GPP-PAR curve (Figure 1, black line).Forcing a GPP model with spatially-averaged PAR could lead to overestimated GPP, since the full spatial distribution of PAR would include values in the light-saturated part of the curve while the spatial average would fall closer to the unsaturated part.In that case, higher values of PAR increase average PAR but do not substantially increase average GPP.
Indeed, we found that GPP calculated from the spatially-averaged PAR was up to 1.5 µmol m −2 s −1 higher (Figure 1, red circle) than the average of GPP calculated from the spatially-variable PAR (Figure 1, blue line).The difference between these two estimates is 24% of the observed standard deviation in daytime half-hourly GPP at a site near Lamont, Oklahoma [43], and indicates the importance of spatial variability in climate drivers of GPP.This likely represents an upper bound on such effects at this site, as the GPP-PAR curve used here has a high degree of nonlinearity to illustrate the utility of our method for a general class of problems.Related problems for which this method is applicable include modeling spatiotemporal variability in soil moisture and leaf area index, and the nonlinear influence of these variables on evapotranspiration and surface energy partitioning [45].

Conclusions
This study presents an easy and practical way of modeling spatial variability in cases where hyper-dimensional covariance (variogram) structures are required to describe observed spatial variability.It could be especially useful in modeling anisotropic conditions with multiple axes of anisotropy, often encountered in environmental and ecological sciences, or in cases where multiple units have to be used to span dimensions of interest (e.g., space and time).The method is simpler than classical approaches to modeling anisotropy (see Chiles and Delfiner, [1], for details) and requires fewer parameters to be estimated from the data.It represents an alternative to approaches that use time series and spatial analysis conjointly to represent variability at unsampled locations (e.g., Romanowicz et al. [46]).
We showed that this approach yields admissible hyper-dimensional covariance models using valid lower-dimensional covariance models as basic blocks.It allows a relatively easy modeling of full 3D-anisotropic space evolving in time, a task that so far has been challenging.
We showed that the original sequential product-sum modeling method can yield results slightly different from the method we propose in this study due to the existence of tolerances and a potentially higher estimated nugget.We simplified the modeling approach and made it applicable to sparse data by proposing an "all at once" modeling approach.
One of the deficiencies of this study is the lack of the specific-case analysis of the criteria for choosing a product-sum model and its higher dimensional derivates instead of other covariance models.While it should be the method of choice in modeling the data characterized by negative nonseparability, it might not be suitable for data characterized by positive non-separability and pointwise non-separability (both positive and negative).For more information on this topic, please see De Iaco et al. [19], De Iaco and Posa [17] and De Iaco et al. [18].
Many ecological processes are nonlinear and can result in biased predictions when using coarseresolution datasets to drive ecological models.We presented a straightforward approach to estimate spatiotemporal variability and demonstrated its application to spatial scaling of ecological processes.This approach relaxes previous assumptions regarding isotropicity of variability for hyper- Due to the saturation of GPP at high PAR, estimates of GPP from the spatially-averaged PAR (red circle) are higher than the actual regional-scale GPP (blue line).

Conclusions
This study presents an easy and practical way of modeling spatial variability in cases where hyper-dimensional covariance (variogram) structures are required to describe observed spatial variability.It could be especially useful in modeling anisotropic conditions with multiple axes of anisotropy, often encountered in environmental and ecological sciences, or in cases where multiple units have to be used to span dimensions of interest (e.g., space and time).The method is simpler than classical approaches to modeling anisotropy (see Chiles and Delfiner, [1], for details) and requires fewer parameters to be estimated from the data.It represents an alternative to approaches that use time series and spatial analysis conjointly to represent variability at unsampled locations (e.g., Romanowicz et al. [46]).
We showed that this approach yields admissible hyper-dimensional covariance models using valid lower-dimensional covariance models as basic blocks.It allows a relatively easy modeling of full 3D-anisotropic space evolving in time, a task that so far has been challenging.
We showed that the original sequential product-sum modeling method can yield results slightly different from the method we propose in this study due to the existence of tolerances and a potentially higher estimated nugget.We simplified the modeling approach and made it applicable to sparse data by proposing an "all at once" modeling approach.
One of the deficiencies of this study is the lack of the specific-case analysis of the criteria for choosing a product-sum model and its higher dimensional derivates instead of other covariance models.While it should be the method of choice in modeling the data characterized by negative non-separability, it might not be suitable for data characterized by positive non-separability and pointwise non-separability (both positive and negative).For more information on this topic, please see De Iaco et al. [19], De Iaco and Posa [17] and De Iaco et al. [18].
Many ecological processes are nonlinear and can result in biased predictions when using coarse-resolution datasets to drive ecological models.We presented a straightforward approach to estimate spatiotemporal variability and demonstrated its application to spatial scaling of ecological processes.This approach relaxes previous assumptions regarding isotropicity of variability for hyper-dimensional analysis, which enables applications across a wider range of environmental and ecological problems.

Software and/or Data Availability Section
The documented Matlab source code (functions) is available at the ResearchGate website [47].The code is made available under 'CC BY' license terms [48].The code was developed by Jovan Tadić in Dec 2017 in Matlab 2016a.Contact details: e-mail: jtadic@lbl.gov.It can be run on a personal computer without special requirements.The archive contains Matlab functions used to model 3D covariance structure using three independent variogram models, presented in the article.The size of the zipped archive is 5 kb.

Figure 1 .
Figure 1.Application of kriged shortwave radiation (i.e., PAR) to estimate gross primary productivity (GPP) over a 2.5° longitude × 2.5° latitude region (~220 km) in central Oklahoma.The black curve illustrates the nonlinear relationship between GPP and PAR from a canopy model.The histogram shows the distribution of PAR values within the region on May 30 of 2006 (16:00-6:30 local time).Due to the saturation of GPP at high PAR, estimates of GPP from the spatially-averaged PAR (red circle) are higher than the actual regional-scale GPP (blue line).

Figure 1 .
Figure 1.Application of kriged shortwave radiation (i.e., PAR) to estimate gross primary productivity (GPP) over a 2.5 • longitude × 2.5 • latitude region (~220 km) in central Oklahoma.The black curve illustrates the nonlinear relationship between GPP and PAR from a canopy model.The histogram shows the distribution of PAR values within the region on May 30 of 2006 (16:00-6:30 local time).Due to the saturation of GPP at high PAR, estimates of GPP from the spatially-averaged PAR (red circle) are higher than the actual regional-scale GPP (blue line).