New parameterized model for GPS water vapor tomography

Water vapor is the basic parameter used to describe atmospheric conditions. It is rarely contained in the atmosphere during the water cycle, but it is the most active element in rapid space–time changes. Measuring and monitoring the distribution and quantity of water vapor is a necessary task. GPS tomography is a powerful means of providing high spatiotemporal resolution of water vapor density. In this paper, a spatial structure model of a humidity field is constructed using voxel nodes, and new parameterizations for acquiring data about water vapor in the troposphere via GPS are proposed based on inverse distance weighted (IDW) interpolation. Unlike the density of water vapor that is constant within a voxel, the density at a certain point is determined by IDW interpolation. This algorithm avoids the use of horizontal constraints to smooth voxels that are not crossed by satellite rays. A prime number decomposition (PND) access order scheme is introduced to minimize correlation between slant wet delay (SWD) observations. Four experimental schemes for GPS tomography are carried out in dry weather from 2 to 8 August 2015 and rainy days from 9 to 15 August 2015. Using 14 days of data from the Hong Kong Satellite Positioning Reference Station Network (SatRef), the results indicate that water vapor density derived from 4-node methods is more robust than that derived from that of 8 nodes or 12 nodes, or that derived from constant refractivity schemes and the new method has better performance under stable weather conditions than unstable weather (e.g., rainy days). The results also indicate that an excessive number of interpolations in each layer reduce accuracy. However, the accuracy of the tomography results is gradually reduced with increases in altitude below 7000 m. Moreover, in the case of altitudes between 7000 m and the upper boundary layer, the accuracy can be improved by a boundary constraint.


Introduction
Based on the GPS meteorology technique (Bevis et al., 1992), two experiments were carried out to assist the development of GPS integrated water vapor (IWV) inversion: GPS/STORM (Rocken et al., 1995) and the GPS-Winter Icing and Storms Project experiment (Gutman et al., 1994).Both experiments showed that GPS is a cost-effective and reliable means of continuously monitoring IWV, with accuracies comparable to water vapor radiometers (WVRs) and radiosondes (RSs).With the development of the GPS-IWV network, many studies have shown that IWV determined by GPS can achieve an accuracy of 1-2 mm (Duan et al., 1996;Emardson et al., 2000;Liang et al., 2015;Niell, 2001;Raja et al., 2008;Rocken et al., 1993).GPS/IWV has been applied to improve the quality of numerical weather prediction (NWP) models and most recently GNSS has been used for nowcasting and help with the predictability of severe weather (De Haan et al., 2009;Gendt et al., 2004;Ichikawa et al., 2012;Rohm et al., 2014;Smith et al., 2007).It is a powerful source of continuous (10-30 min temporal resolution) integrated water vapor determination for climate studies, weather prediction, and hazard mitigation.However, IWV is a measure of the total amount of water vapor above a certain station, and it cannot provide information on the vertical distribution of water vapor.
In order to meet the demand, GNSS (Global Navigation Satellite System) water vapor tomography has appeared as a promising method of providing information on the fourdimensional distribution of the water vapor in the troposphere.By using extensive slant wet delay (SWD) data collected from a network of GNSS receivers, a four-dimensional humidity field with high spatial and temporal resolution can be reconstructed from tomography.Over the past decade, numerous studies demonstrated the potential of GNSS tomography to retrieve 4-D humidity fields.Unfortunately, because N. Ding et al.: New parameterized model for GPS water vapor tomography the quality of the reconstructed profiles is limited by the constellation of GNSS satellites, the geographic distribution of ground-based receivers, and observation errors (Chen and Liu, 2014;Rohm et al., 2014;Shangguan et al., 2013;Yao et al., 2016), some voxels may not be crossed by any signal during a tomographic process.Consequently, this will lead to an ill-posed inverse problem with incomplete input data.
In this paper, we introduce a parameterized approach based on horizontal IDW interpolation, in which only the vertical constraint is used to ensure that the model is compliant with the features of water distribution in the troposphere.The new algorithm is analyzed using data from SatRef.The experiment mainly analyzes and discusses the influence that different numbers of voxel nodes have on the results of GPS water vapor tomography.Moreover, we address situations in which some voxels are not crossed by any signals.However, grid points not used in any interpolation should not generally be avoided.In fact, this case often occurs in lowerlevel layers with the four-node method, and occurs less often when many points are included in the interpolation.We also address "empty" grids using the inverse distance weighted (IDW) interpolation.The values of an "empty" grid are estimated by calculating the "non-null" grids around it.

Tomography observations
Atmospheric refraction delay, including ionospheric delay and neutral atmospheric delay, is a major source of GPS positioning errors.Ionospheric delay (dispersive delay) can be corrected by using dual-frequency linear combinations.Neutral atmospheric delay, independent of frequency, is an area of great research interest in geosciences and meteorology.
The slant total delay (STD) can be expressed by (Bevis et al., 1992 In the first term, N is the atmospheric refractivity and s is the path of the GPS signal propagation in the troposphere.The second term, [S-G], is due to ray bending.Fortunately, the effect of this geometric delay, which has only a minor effect on STD, can be neglected (Ichikawa et al., 1995).
The delay caused by atmospheric refractivity contains the slowly varying delays caused by dry air and the rapidly varying delays caused by water vapor.This characteristic can be represented by the following formula (Thayer and Gordon, 1974): where P d is the partial pressure of dry air (mbar), P w is the water vapor partial pressure (mbar), T is the temperature (K), and Z −1 d and Z −1 w are the inverse compressibility factors for dry air and water vapor, respectively.In our study, both factors are assumed to be constant (Z −1 d = Z −1 w = 1).The empirical coefficients t 1 = 77.6,t 2 = 64.8, and t 3 = 3.776 × 10 5 K 2 mbar −1 refer to (Bevis et al., 1994).SHD is the hydrostatic portion of the delay in the slant direction.It can be calculated by applying a modified version of Saastamoinen's formula (Saastamoinen, 1972).SWD is the slant wet delay; it can be written as a function of azimuth (φ) and elevation (e): where ZWD is the zenith wet delay, which can be obtained from zenith total delay (ZTD) minus zenith hydrostatic delay (ZHD).m w (e) is the wet mapping function, and G w N and G w E are the wet delay gradient parameters in the north-south and east-west directions, respectively.R is the residual unmodeled delay between the GPS satellite and the receiver.The "zero differences" (ZDs) method (Alber et al., 2000) has been proposed for determining slant-path delays from double-difference post-fit residuals.However, this method has inherent limitations for determining non-homogeneous atmospheres, because the method spreads the inhomogeneous signals over all parameters estimated in the leastsquares step (Elosegui and Davis, 2003).Thus, in our study, R will be removed from the SWD calculation.3 Modeling the humidity field with IDW parameterization IDW parameterization assumes that the troposphere is divided into a number of vertical layers (see Fig. 1).SWD produced by one signal (blue line) is discretely modeled as the sum of the values at the points of intersection between straight lines and vertical layers, multiplied by the corresponding distance.There are no vertical variations assumed within a layer.In our study, the true path of the signal between the satellite and the receiver will be replaced by a straight line, and the vertical structure of nonuniform layers is used.
The SWD of a signal propagation path can be expressed by the formula: where n is the number of vertical layers, P k (x,y,z) is the value of the wet refractivity (mm km −1 ) at position (x, y, z) of the intersection at the kth layer, and l i is the distance traveled by the ray in the layer.
The IDW parametrical method explicitly implements the assumption that the values of wet refractivity in close proximity are more alike compared to those that are farther away.It is clear that the wet refractivity at an intersection can be calculated based on a weighted average of the values at voxel nodes that surround it.The most common method to define intersection value is to utilize a limited number of voxel nodes on the same layer to express it.For example, in Fig. 2, four adjacent nodes (x 1 node , x 2 node , x 3 node , x 4 node ) in close proximity to intersection P will be used to estimate the value of wet refractivity at point P above the reference stations.
Based on IDW interpolation, P k (x,y,z) can be expressed by the following function: where D(P , x i node ) is the distance between pierce point P k (x,y,z) and voxel node x i node , x i node is the value of water vapor at a certain voxel node (including the four adjacent nodes close to point P in the case depicted in Fig. 2), i is the node number that will be used for IDW interpolation, and m is the total number of x i node used for interpolation.Therefore, SWD i (generated by signal propagation in the troposphere) can be given as where n is the number of vertical layers, k is the sequence number of each layer, and f (x i node ) is the function of x i node introduced in Eq. ( 8).Equation ( 7) can be written in vector form.Combining all SWD i vectors together, Eq. ( 7) becomes where vector SWD represents the SWD observations vector, x node is the state vector of the wet refractivity at all design voxel nodes, and C is a mapping matrix, which defines the mapping of x on the SWD observations.In order to reconstruct the characteristics of water vapor distribution in a vertical direction, the vertical constraint will be utilized.We will assume an exponential law (Davis et al., 1993), which can be used to represent an average water vapor profile: where n k is the wet refractivity of the kth layer, h k is the height of the kth mlayer, and H is the water vapor scale height, which can be calculated using Eq. ( 9) (Zhang et al., 2015): where W is the total vertical water vapor content in g m −2 and ρ s is the surface humidity in g m −3 .W can be obtained from PWV.To support meteorological monitoring, groundbased GPS networks provide meteorological parameters that can be used for calculating ρ s .Based on Eq. ( 10), the vertical constraint will be constructed by where V is the vertical constraint coefficient matrix; x node is the same as that occurring in Eq. ( 6).By combining all the SWD observations (Eq.7) and vertical constraints (Eq.10), the solution model of IDW tomography in a matrix form is given as where A is the coefficient matrix of an IDW tomography model; Eq. ( 11) can be written as Ax = b.
4 Tomographic reconstruction using ARTs based on PND access order scheme Algebraic reconstruction techniques (ARTs) have been successfully used to reconstruct the humidity field (Bender et al., 2011;Chen et al., 2014).An advantage of ARTs is that they have high numerical stability even under adverse conditions.Moreover, it is relatively easy to incorporate prior knowledge into the reconstruction process.
The solution model of IDW tomography is defined by Ax = b (Eq.14).It is common to compute an approximation of the tomographic solution using an ART via the following formula: where a i is the ith row of A and b i is the ith row of b, x k is the kth iterative solution (which will be used to update the next solution x k+1 ), and λ is the relaxation parameter; its empirical value is 0.2 in our study.The classic method of the ART family is Kaczmarz's algorithm (Kaczmarz, 1937); its k iteration consists of a traversal through the m rows of A in ascending order (i.e., from 1 to m in Eq. 11).
There are other ARTs, such as the symmetric method (Björck and Elfving, 1979) and the randomized method (Strohmer and Vershynin, 2007), that can distinguish according to the order in which the rows are processed.Therefore, the access order of the ART has a significant effect on its practical performance (Herman, 1980).The importance of access order has been recognized in medical applications of ART (Mueller et al., 1997;Herman and Meyer, 1993;Hounsfield, 1976;Robb et al., 1974).
In our study, an access order scheme based on prime number decomposition (PND) will be proposed.It is desirable to order the SWD observations such that subsequently applied SWDs are largely uncorrelated (Ding et al., 2016).This means that consecutively available SWDs must have significantly different values, because the value of SWD is determined by the azimuth and elevation angles of a signal.If the SWDs in a set have similar values solved by ART, the results tend to move away from the desired solution, which delays convergence.To summarize, the principle of SWD access order is that in a subsequence of iterations of Eq. ( 15), steps should be as independent as possible from the previous steps.
Following the decorrelation principle, the PND access order scheme (Herman and Meyer, 1993) is presented in this section.In the first step, we sort all elements of the SWD vector (in Eq. 10) in descending order; the sorted equations will be numerated from 0 to M. The mapping relationships between the PND access order of equations and prime factors m can be established by the PND access order scheme.As a simple example, we set M = 473 (M + 1 must not be a prime number).The detailed steps to produce the PND access order in this case are given by Table 1.Because of space limitations, only the top 30 PND permuted sequences of equations are displayed.It is clear that if we subdivide the permuted sequence into pairs, there are 237 between them (the largest possible value).If the permuted sequence is subdivided into groups of six, the common difference of 79 arithmetic progressions is obtained, which is also as large as possible.This shows that the PND access order must comply with the decorrelation principle.5 Water vapor tomography from the Hong Kong network

Network and voxel node model
SatRef is a local satellite positioning system covering the extent of Hong Kong.The network consists of 18 continuously operating reference stations (CORSs); 14 of them were used in this study.The horizontal and vertical station distributions are presented in Fig. 3a and b, respectively.The area of investigation ranges from 113.749 to 114.474 • E longitude and from 22.115 to 22.651 • N latitude (Fig. 3a) with a height range of 0-10 800 m (Fig. 3b).The number of spatial voxel nodes (Fig. 3a) used for tomography is 6 × 4 × 10 with a horizontal resolution of 0.145 • (15 km).Nonuniform layers are used in the interval from 500 to 3800 m (Fig. 3b).We hypothesized that the value of water vapor at any point in the layer is determined by a weighted average of the value at the voxel nodes closest to the point.
Based on the number of voxel nodes x i node (Fig. 2) used for IDW interpolation points P (Fig. 2) and the classic method (Flores et al., 2000;Hirahara, 2000) to make refractivity constant within a voxel, four solutions are presented: -4nodes: four voxel nodes (neighbors) from the voxel node model are used when calculating a wet refractivity value for the piercing point.
-8nodes: eight voxel nodes are used when calculating a wet refractivity value for the piercing point.-12nodes: 12 voxel nodes are used when calculating a wet refractivity value for the piercing point.
-Constant: the refractivity is constant within a voxel when calculating a wet refractivity for the piercing point.

Analysis of total statistics
The four solutions are compared with radiosonde (RS) data recorded in Hong Kong.Because of limited space, only a representative example from 9 August 2015 (DOY 221) 00:00 and 12:00 UTC is shown (see Fig. 4).It is clear that all of the water vapor profiles (red lines) reconstructed by the tomographic solutions accord with the radiosonde data.At 00:00 UTC, all profiles and the scatter of radiosonde data decrease exponentially with increases in height.However, at 12:00 UTC, radiosonde data indicate disturbances in water vapor density and are thus less stable than at 00:00 UTC.In each layer (Fig. 3b), the mean water vapor density is computed.It indicates that it is difficult to retrieve humidity data with high vertical resolution (e.g., radiosonde) because of the limited accuracy of the tomography solution.However, consistency is maintained with the changing trends of the radiosonde data.It is difficult to determine which solution is best suited for GPS water vapor tomography from the statistics of these data (Table 2).These solutions have their respective advantages.The 4nodes parameterization presents more accurate results than others in terms of RMSE.The interquartile range (IQR), to some extent, reflects the discreteness of datasets.Based on the IQR values (Table 2) and the water vapor density scatter plots (Fig. 5a-d), it is clear that the tomographic solutions computed by 4nodes parameterizations are more concentrated than those computed by the others.However, compared with the other methods, the 12nodes tomography has the smallest bias (Table 2) and the highest discretization.
The statistical characteristics of the differences between the four schemes and the radiosonde data are also presented by box plots (Fig. 5e).They contain five characteristic values: the first and third quartiles (Q1 and Q3) are located at the bottom and top of the box; the second quartile (Q2) is located inside the box (the median) and at the ends of the whiskers (up- per and lower bound) at Q1 − 1.5(IQR) and Q3 + 1.5(IQR).In addition, whiskers are applied to detect outliers (cross in Fig. 5e).All of these are summarized in Table 3.
IQR is the difference between the Q3 and Q1 quartiles.It is important because IQR represents the spread of data and, unlike the total range, it is not affected by outliers.Q2 is the measure of central tendency and is in good accordance with the bias (as the statistics in Table 2 indicate).Upper and lower bounds can be used to identify outliers.The constant scheme has the smallest number of outliers, but it also has maximum bounds.The 8nodes and 12nodes schemes have similar bounds.However, if we use the bounds of 4nodes, the numbers of outliers in the 8nodes, 12nodes, and constant schemes are 66, 30, and 30, respectively.Therefore, comparison of 4nodes parameterization with the other schemes shows that it has a relatively minimal number of outliers.

Analysis of statistics and relative error in each layer
The above analyses are based on total statistics of data from a 7-day period.However, the precision of water vapor tomography is highly influenced by the vertical structure of the humidity field.For example, the same tomographic error value in different layers has different qualities.Therefore, total statistics cannot provide a comprehensive analysis of the accuracy of tomographic reconstruction.The relative index is introduced for tomographic stratification analysis in Fig. 7.We adopt nonuniform layers within the interval from 500 to 3800 m.The box plots (Fig. 6a-d) reveal that the errors of 4nodes and 8nodes in each layer are, in general, more concentrated than those of 12nodes and the constant method.
The results accord with the total statistical analysis, which was illustrated in a previous section.A comparison of the three parameterizations in terms of RMSE (Fig. 6e-h) and relative error (Fig. 6i-l) shows that 4nodes has relatively small values for these two parameters.However, by observing the changes between RMSE and relative error, opposite tendencies emerge in total parameterizations.For example, between 4000 and 7000 m, the RMSE of the constant method decreased from 1.749 to 1.078 g m −3 , a decrease of approximately 38 %.However, the average humidity between 4000 and 7000 m declined by 77 %, which led to a relative error increase from 39 to 62 %.This illustrates the following: with increases in height, the decrease in moisture content makes it difficult for water vapor tomography techniques to retrieve humidity in high altitudes.However, the changes in 4nodes' values caused by increases in height were smaller than those of the others.The voxel nodes in the top layers are assigned a fixed constant value (0 g m −3 ).This is a widely used water vapor tomography method (Flores et al., 2000) that markedly reduces the relative errors in the upper boundary layer.

GPS water vapor tomography in dry weather
Four methods are used again in this period.All of the water vapor profiles reconstructed by the tomographic solutions are in good agreement with the radiosonde data.To avoid repetition of conclusions, this section focuses on the difference between experiments during rainy days  and dry weather .The results of RMSE in experiment 2 (Table 4) increase accuracy by more than 2 times compared to that in experiment 1 (Table 2).In terms of bias and IQR, consistency is maintained with the case of RMSE.This means that the new method has better performance under stable weather conditions than unstable weather (e.g., rainy days).However, by comparing the results of outliers in both experiments, the main reasons for outliers in Table 4 having a remarkable increase are due to experiment 2 having smaller bound than those of experiment 1.
Comparing four methods in experiment 2, 4nodes has the smallest value of RMSE and number of outliers.Similarly, 4nodes is more concentrated than other methods.This case is reflected in the value of IQR (|Q1-Q3|) and Q2 (Table 4).Interestingly, the 12nodes tomography has the smallest bias again with the highest discretization.Outliers, which are the values out of the bounds (upper and lower bounds in Table 4), allow, to some extent, for judgement regarding the stability of tomographic methods, and 4nodes is better than other methods.RMSE and relative error in each layer (Fig. 7) have good agreement with experiment 1.It is show that as altitude increased, opposite tendencies emerged between the RMSE of each layer and the corresponding relative error.In the top layers, a boundary constraint condition dramatically reduces the RMSE and relative errors.

Conclusion and outlook
In this paper, a new GPS tomographic parameterization approach based on IDW interpolation is proposed.This approach can reconstruct a humidity field without using horizontal constraints and prevent situations in which some voxels are not crossed by satellite paths.On the other hand, instead of dividing the troposphere into several layers with identical heights, the vertical structure of the tomography model adopts nonuniform layers to satisfy inherent characteristics of water vapor distribution and to lower the effect of the difference magnitude between the calculated tomographic results.In order to minimize correlation in projection access, a PND access order scheme is developed to order the SWV observations such that subsequently applied SWV values are largely uncorrelated.Based on the number of voxel nodes selected for IDW interpolation, and the constant method that makes refractivity constant within a voxel, four schemes are designed to retrieve water vapor density into voxel nodes.In addition, a vertical constraint is adopted to ensure the characteristics of vertical water vapor distribution.
Tomographic experiments using GPS data collected over Hong Kong from DOY 214 to 227 2015 validate our proposed GPS tomography-based approach.We discuss and analyze 4nodes, 8nodes, and 12nodes methods as well as a constant method.For the overall dataset, the 4nodes method offers the highest accuracy compared to the other three meth- ods, with perturbations of 1.627 g m −3 on rainy days and 0.826 g m −3 in dry weather, respectively.Moreover, the new parameterized method has better performance under the stable weather conditions than the unstable weather (e.g., rainy days).However, owing to the limitations of tomographic accuracy, it is more difficult to use the mean water vapor density computed by GPS tomography to retrieve the humidity field, compared with using the high vertical resolution capability of the radiosonde.Based on the results of different statistical methods, we can conclude that 4nodes parameterization has a relatively small number of outliers, and that the errors it causes are more concentrated than in the other schemes.This indicates that tomographic results derived from 4nodes offer higher stability and reliability.Tomographic layering is studied, and results show that as altitude increases, the opposite tendency between the RMSE of each layer and the corresponding relative error can be observed.It indicates that because the moisture content decreased exponentially with increases in altitude, it was increasingly difficult to retrieve humidity using water vapor tomography.
However, in the upper boundary layer, relative error is reduced markedly with the use of boundary constraints.
In future studies, the horizontal structure of the humidity field must be improved by adjusting the node position in each layer to fit the distribution of satellite rays.Flexible layout is the advantage of the voxel node model.In the future, the fusion of GNSS and external measurements from other sensors in the GPS tomography system will be a potential means to enhance the stability and reliability of water vapor tomography and to decrease tomography intervals.

Figure 1 .
Figure 1.Vertical structure of IDW tomography model and the signal rays (blue line) crossing the vertical layers.The value of SWV divides into the value of water vapor at pierce points (red points).

Figure 2 .
Figure 2. A case showing an IDW tomography model in a certain layer.Black points are voxel nodes used for predicting the value of the pierce point (red point); the gray wireframe indicates the layer used in the model, and the length of the double arrows represents the distance between the point P and voxel node x node .

Figure 3 .
Figure 3. (a) Distribution of the Hong Kong reference stations (red triangle) and the King's Park Meteorological Station (red square).The gray wireframe indicates the layer that contains voxel nodes (black dot) used in the tomographic processing.The study area was discretized into 6 × 4 × 10 voxel nodes for the water vapor tomography.(b) Vertical structure of the voxel node model used in the tomographic processing and height distribution of the GPS stations (black spots).
Two experiments are carried out in this paper.The first experiment (experiment 1) is used for reconstruction of the humidity field on rainy days from 9 to 15 August 2015 (day of year (DOY) 221-227).Thunderstorms continued to affect the experimental region (Hong Kong) during this period, and no rain occurred during the remainder of August 2015.During the rainy days, the moisture content increased and changed dramatically in the troposphere; this was suitable for verifying the feasibility of our new GPS water vapor tomography method.For completeness another experiment (experiment 2) with the same number of days is conducted in dry weather from 2 to 8 August 2015.With plenty of sunshine, conditions became very hot during the day, with maximum temperatures exceeding 33 • C on 3-7 August.The summer heat grew even more intense on 8 August.Temperatures at the observatory soared to a maximum of 36.3 • C on the afternoon of 8 August, an all-time high since records began in 1884.In the case of the doubledifference process model, the absolute values of atmospheric parameters must be obtained by introducing baselines longer than 500 km(Rocken et al., 1993).Three IGS sites (BJFS, PIMO, and LHAZ) were added to the Hong Kong network as auxiliary stations to meet the long baseline requirement.Moreover, radiosonde data from King's Park Meteorological Station (HKKP) are used to assess the quality of GPS water vapor tomography.The radiosonde collects water vapor data with high a vertical resolution, and the sounding balloons are launched twice a day at 00:00 and 12:00 UTC.The available radiosonde data from HKKP are obtained for comparison with the results of GPS water vapor tomography.The wet refractivity estimated by solving tomography equations will be further converted to water vapor density.The relation-

Figure 5 .Figure 6 .
Figure 5. Tomographic solutions compared with radiosonde data.(a-d) Scatter plots of water vapor density with 4nodes, 8nodes, 12nodes, and constant parameterization and (d) box plots of the differences between the four methods and radiosonde data at 00:00 and 12:00 UTC from 9 to 15 August 2015.

Figure 7 .
Figure 7. Statistics of the humidity differences between tomographic solutions and radiosonde data in each layer from DOY 221 to 227 2015.(a-d) RMSE and (e-h) relative error of the humidity difference in 10 layers using 4nodes, 8nodes, 12nodes, and constant parameterizations.

Table 1 .
Prime number decomposition access order.

Table 2 .
Statistics showing the differences between the four schemes and water vapor measurements constructed by HKKP's radiosonde data for 7 daysin 2015.

Table 3 .
Characteristic values of box plots of the four schemes and water vapor measurements constructed by HKKP's radiosonde data for 7 days (DOY 221-227) in 2015.

Table 4 .
Statistics and characteristic values of box plots showing the differences between the four schemes and water vapor measurements constructed by HKKP's radiosonde data for 7 days (DOY 214-220) in 2015.