Modelling quasi-periodic signals in geodetic time-series using Gaussian processes


 Seasonal signals in geodetic time-series have long been recognized to be associated with environmental phenomena such as polar motion, atmospheric loading, groundwater loading and other hydrological processes. Modelling these periodic signals is crucial for the geophysical interpretation of these time-series. The most common approach used for resolving seasonal (annual and semi-annual) signals is their approximation by sinusoidal functions with constant amplitudes. However, because of their environmental source, seasonal signals are likely to be quasi-periodic. In this study, we investigate a Gaussian process (GP) to model quasi-periodic signals in geodetic time-series, a flexible method that allows capturing the variability structure in the data using covariance functions. We use the Markov Chain Monte Carlo method to evaluate the posterior probability density function. To test its effectiveness, we apply this method to a synthetic time-series in the presence of time-correlated noise. We find that the GP model provides a better fit to the time-series, resulting in time-series residuals with fewer systematic effects. We use the GP model to estimate the secular velocity of selected GPS sites from Antarctica and Alaska, as well as an example of Gravity Recovery and Climate Experiment time-series. The Bayesian aspect of the GP model allows inferring the linear velocity ensemble in the vicinity of the true solution while taking into account the quasi-periodic systematics in the time-series.


I N T RO D U C T I O N
Modern space geodetic techniques offer surface position time-series with an unprecedented level of precision. Coordinate time-series derived from the Global Positioning System (GPS) are routinely used to define terrestrial reference frames and to study of a variety of geophysical phenomena, including crustal movements related to plate tectonics (e.g. Feigl et al. 1993), glacial isostatic adjustment (e.g. Milne et al. 2001), magmatic processes (e.g. Segall 2010), atmospheric surface loading effects (e.g. Van Dam et al. 2001) and hydrological surface changes (e.g. Borsa et al. 2014). Surface displacement time-series contain three main types of signals: (i) long-term slowly changing signals such as linear secular motion associated with the interseismic motion of plate tectonics; (ii) transient changes such as sudden position jumps related to co-seismic displacement or post-seismic deformation and (iii) periodic oscillations corresponding to the crustal response to seasonal environmental changes. Although the first two types of signals were identified since the early days of the GPS development, the periodic oscillations (annual and seasonal) have only been adequately studied in the last two decades after the deployment of multiple permanent GPS stations offering continuous records of the positions, which was not possible in the case of campaign GPS measurements. The annual and semi-annual signals are the dominant observed periodicities in position time-series and they are often assigned to local and/or global environmental effects, such as thermoelastic effects, atmospheric loading and groundwater pumping. Systematic positioning errors, as well as reference frame mismodeling, can also contribute to the position signals (e.g. Ray et al. 2008;King & Watson 2010). The impact of the annual signal on the estimation of secular velocities from geodetic time-series is today well known since it was first demonstrated by Blewitt & Lavallee (2002). If not accounted for, it can induce a significant bias in the estimation of the linear velocity. Therefore, most geodetic studies today use models that include sinusoids with annual frequency terms to reproduce the annual/seasonal cycles. Although this approach might seem sufficient for capturing the annual periodicity in GPS time-series, it suffers from a limitation resulting from the assumption that the amplitude of the signal is constant and cannot change from year to year. As the main sources of the annual and seasonal oscillations in GPS timeseries are environmental related, it is expected that they exhibit non-stationary behaviour. Therefore, a more complete description of the time-series model would include deviations from the purely repeating annual cycle, which we refer to here as quasi-periodicity.
A number of studies have investigated the quasi-periodic variability when modelling GPS time-series analysis. Langbein (2004) proposed a bandpass filter to model noise sources where the 1 cycle yr −1 frequency leaks into adjacent frequencies. The passband for the filter considered lies between 0.5 and 2 cycles yr −1 . The sharpness of the filter was adjusted by the number of poles used to construct the filter. He showed that approximately 30 per cent of the electronic distance meter (EDM) baseline seasonal noise spread over a frequency band around 1 cycle yr −1 . Bennett (2008) proposed a semi-parametric method to investigate the effect of quasi-periodic signals on secular velocity. He showed that if the temporal variations of the amplitude of the annual motion are not accounted for, the residual-error of time-series may contain significant power at the semi-annual band. Davis et al. (2012) used a stochastic approach using a Kalman filter to model the seasonal variability. A more general version of the Kalman filter has been proposed by Didova et al. (2016), where the noise parameters were estimated using a least-squares approach. Chen et al. (2013) applied the singular spectrum analysis method to model time-varying signals in weekly GPS time-series. Klos et al. (2018) published a complete summary and a comparison between published methods of modelling time-varying seasonal signals in GPS time-series. Using synthetic time-series, they showed that the precision of the recovery of quasi-periodicity depends on the annual-signal-to-noise ratio. They also showed that the Kalman Filter approach outperforms other published techniques by modelling between 77 and 90 per cent of the total variance of the varying seasonal signal. However, this can only be achieved by a special tuning of the noise parameters because of the non-convex nature of the likelihood optimization problem.
Given the importance of the noise characterization in the modelling and interpretation of the time-varying signals in GPS timeseries, a more sophisticated method would involve a simultaneous estimation of the trajectory model of the surface displacement of the GPS site along with systematic artefacts and the noise processes. In this paper, we explore a different approach using the Gaussian process (GP) to model quasi-periodic signals in GPS coordinates timeseries. GPs are widely used in machine learning and have been used in different aspects of signal processing (Rasmussen & Williams 2006). GPs are Bayesian non-parametric methods used for regression and classification problems. Despite the fast-growing interest in using GPs in many physical applications, very few attempts have been made to use them for geophysical problems. Hines & Hetland (2018) used Gaussian process regression (GPR) to detect transient strain resulting from slow-slip events. Sarkar et al. (2019) used GPs to predict fast tidal currents. Valentine & Sambridge (2020) described a theoretical framework of using GPs in geophysical inverse problems. Here, we use GPR to infer the secular rate of GPS time-series in the presence of quasi-periodic signals and correlated noise. We begin by introducing the GP model and the Bayesian estimation of the parameters. In Section 2, we apply our GP model using synthetic time-series. Then we show an application to real time-series.

Linear trajectory model
The time evolution of a GPS station coordinates y(t) can be described as a trajectory model in the form of a linear combination of a secular linear velocity and a seasonal signal: where t R is a reference time, y R is a reference position, v is a timeconstant velocity, N f is the number of frequencies, s k and c k are the Fourier coefficients and ω k is the angular frequency. In terms of a matrix formulation, the above eq. (1) can be re-written as where the coordinate observations are stored in vector y, β represents the model parameters and ζ represents the residuals assumed to be independent Gaussian noise. The X matrix is an n × m matrix where n is the number of observations (epochs) and m is the number of model parameters. This model is linear in its coefficients, and so it can be solved using a linear least-squares method:

Gaussian process model
Instead of considering a deterministic form of the trajectory model, as described above, we model the coordinate time-series as a GP, where the joint distribution probability of the observations y is a multivariate Gaussian distributed about a mean function f with a covariance matrix , given by where φ are the parameters of the mean function and θ are the covariance matrix parameters referred to as hyperparameters. The covariance matrix contains information about the observation's systematics and the noise. For example, white noise would be represented by a matrix with variances in each element of the diagonal. The elements of the covariance matrix are formed in the GP problem using covariance kernel functions k(t i , t j ): The term σ 2 i δ i j represents the white noise, and δ ij is the Kronecker delta. In our case, σ 2 i are taken to be the uncertainties from the time-series position solution.
The main advantage of the GP method is its flexibility in modelling complex signals which may include quadratics, sinusoids or complex shapes without prior knowledge of the signal form. To do this, a parametrization of the covariance matrix is needed. There are various forms of the kernel functions in the literature (Roberts et al. 2012) and their choice usually arises from the physical nature of the problem treated. (Rasmussen & Williams 2006) proposed the product of a periodic function and a squared exponential function to model the decay from exact periodicity of the Mauna Loa CO 2 concentration. The Matérn-3/2 family functions is an alternative kernel flexible enough to model rough local variations without the addition of hyperparameters. We opted for the following kernel (Foreman et al. 2017): (6) where σ 2 and ρ are the amplitude and the timescale of the variations τ = t i − t j . In the limit → 0 this covariance function becomes the Matérn-3/2 function (Rasmussen & Williams 2006). The parameter controls the quality of the approximation. As indicated by Foreman et al. (2017), the value of epsilon will depend on the specific data set and the precision requirements for the inference. We have Downloaded from https://academic.oup.com/gji/article/226/3/1705/6253216 by guest on 30 July 2021 tested different values of epsilon between 0.0001 and 0.001 using the same synthetic time-series described later in Section 2 and we obtained similar results. We tested the stochastically driven damped harmonic oscillator kernel (SHO; Foreman et al. 2017); however, our Markov Chain Monte Carlo (MCMC) sampler did not converge. The choice of Matérn 3/2 kernel was motivated by its flexibility in allowing the description of short-term temporal variations as well as rough aspects in the coordinate signals, while using only two hyperparameters. Additionally, we model the mean function,f(t), as a linear trend for the secular long-term motion, and a sinusoidal function, which represents the constant part of the annual periodic signal: where d is the y-intercept, v is the secular velocity, and A and B are the amplitudes of the annual sine and cosine functions. The Heaviside step function H is included only if an offset (point-change) due to an instrument change or a co-seismic displacement is identified at an epoch t j (n j is the total number of offsets). The log-marginal likelihood of the GP is written as where N is the number of data points and r is the residual vector from the mean function The calculation of the posterior log-likelihood (eq. 8) requires the evaluation of the inverse and the determinant of matrix . This makes the computational cost proportional to N 3 . Therefore, the evaluation of the GP likelihood using MCMC methods becomes very expensive, especially for large time-series (more than 10 4 points). A remedy of this limitation was proposed by Ambikasaran et al. (2016), where the matrix can be hierarchically factored into a product of block low-rank updates of the identity matrix. The application of this method is restricted to one-dimensional datasets and the family of available kernels is also limited, although a certain number of popular kernels can be approximated with the combination of CELERITE kernels (Tyler et al. 2020). In our study we used the CELERITE 1 python package for the implementation of our GP. The hyperparameters θ , including the mean function parameters φ, were obtained using a full Bayesian inference via an MCMC method. The sampling of the joint posterior probability was performed using the emcee 2 package (Foreman-Mackey et al. 2013). We ran 500 steps burn-in and 2000 steps of MCMC. We limited the explored space of the hyperparameters based on our physical understanding of the position signals to prevent the MCMC chain from moving into unrealistic parameter space. The Bayesian approach allows us to infer the ensemble of optimal models that can explain observed data.

T E S T O N S Y N T H E T I C T I M E -S E R I E S
We used a synthetic position time-series to demonstrate the performance of the GP method proposed in the previous section. We used the same kinematic model as in Bennett (2008), including a long-term linear rate (5 mm yr −1 ) in addition to a constant annual signal and a sinusoidal variation representing the time-variable part of the quasi-periodic signal: where A and B represent the time-averaged signal amplitudes (2.4 mm) and c(t) is the deviation function from the time-averaged signal. Here we define c(t) as a sinusoidal function with 5 yr period and an amplitude similar to the time-averaged. We consider the time-variable part of the signal in phase with the time-averaged part, such that p = tan −1 (b/a). An offset of 3 mm was added at the 3.5 yr epoch to simulate a typical instrumental change in GPS time-series. We added a noise consisting of a combination of white (1 mm) and flicker (1 mm yr −1/4 ) models. We chose a time-series with 7 yr, well beyond the minimum data span recommended by Blewitt & Lavallee (2002) to avoid the degradation in precision arising from the estimation of extra annual signal parameters.To compare the GP model fit with the traditional approach, we fitted the same synthetic time-series using a time-constant seasonal signal. The noise was parametrized using a combination of Generalized Gauss Markov (GGM) and white noise models and the parameters were estimated using maximum likelihood estimation (MLE) in the Hector software (Bos et al. 2013). This will be referred to as the standard model in this study. Additionally, we included in the comparison two other methods that take into account the time-varying seasonal signal: (i) The Wiener Filter-based (WF) method proposed by Klos et al. (2019) and (ii) the Basis Functions (BF) method described in Bennett (2008). We used the Hector implementation of both methods. The WF method consists of using a first-order autoregressive process to model the time-variable seasonal signal, then applying a Wiener Filter adapted to the noise level of data. The BF approach uses orthogonal basis functions to describe the slow varying seasonal signals. Here we use the Hector implementation of this method using the Chebyshev polynomials of a degree 3. Fig. 1 shows the best fit using the three approaches as well as the posterior mean of the GP model. Based on a visual inspection, it is clear that the methods considering the time-varying seasonal signal provide a better fit to the data than the standard approach. The residuals from the standard model show smooth periodicity as a result of not accounting for the deviations from a constant annual signal. Using the BF method produces smaller residuals (Fig. 1d), but with a clear pattern of seasonality. Although the WF method shows a better performance (Fig. 1e), the obtained residuals show a long-period signal indicating that the seasonal periodicity was not fully modelled, unlike the GP model that provides residuals with no periodic systematics (Fig. 1f). The covariance matrix parametrization of the GP model gives considerable flexibility allowing the GP model to represent the quasi-periodic signals in the time-series. Fig. 1(b) shows that the inadequate modelling of the annual signal, in the classic model, biases the estimation of the offset amplitude, similar to the finding of Montillet et al. (2015) who showed a bias of up to 0.5 mm in the co-seismic offsets for an amplitude of the seasonal signal greater than 1 mm. Fig. 2 is the corner plot of the marginal posterior distributions of the hyperparameters. It shows that the distributions are consistent with the true values used in the simulation with a slight correlation between the linear rate and the intercept parameters. Fig. 3 shows the power spectral density (PSD), computed using the Welch method, for the residuals of each model. At the high-frequency level, the power spectra are relatively similar; however, for the frequency band between 1 and 6 cpy, the classical model shows more power than the GP model (Fig. 3a, b) because the standard model ignores the time-dependent annual variations. The PSD of the standard model fits well a power law plus white noise model with a spectral index of −0.99, in line with what has been previously shown by different studies of the noise properties in GPS time-series (Mao et al. 1999;Williams 2003;Langbein 2004;Bos et al. 2008). However, the PSD of the GP model residuals is more flattened with a lower spectral index (κ = 0.2), indicating that the GP models absorb some of the noise in addition to the stochastic annual variability.
To understand the effects of the amplitude variations of the annual signal on the estimation of the secular rates, we compared the velocity estimates from the GP method against the standard model described above, using different time durations of the time-series. We generated 200 synthetic time-series with durations more than 2.5 yr at integer-plus-half year spans to ensure a minimal bias velocity due to constant annual variations. We created 200 time-series for each period of 3.5, 5.5 and 10.5 yr, consisting of a noise model similar to the example described above. Table 1 summarises the results of the comparison. The overall comparison shows that the Median Absolute Difference (MAD) with respect to the true rate decreases as the time span of the time-series increases. For short time spans (3.5 and 5.5 yr), the BF method produces the largest MAD even larger than the standard model with a time constant seasonal signal model. Klos et al. (2018) found similar result when comparing with weighted least-squares for synthetic time-series with low noise amplitudes. The WF performed better than the BF and the standard methods for short time-series, but resulted in comparable misfit for longer time-series. As mentioned by Klos et al. (2019), the misfit reduction also depends on the amplitude of simulated noise. They showed that the WF outperforms when high amounts of noise exist in time-series. The GP method provided the smallest MAD for the three time-series periods, but resulted in the smallest trend uncertainties compared to the other methods. This is possibly related to the fact that the GP method absorbs too much power for periods larger than 10 cpy (Fig. 3). We attribute this to the choice of the covariance. Although the Matérn 3/2 is a flexible function with few hyperparameters, other combinations of covariance functions might provide better modelling with more realistic velocity uncertainties. We leave this question for future investigations.

Antarctic GPS
Different studies have attempted to model the annual signal from GPS time-series by modelling different contributions from loading effects and/or from the surface elastic deformation derived from the Gravity Recovery and Climate Experiment (GRACE; Jiang et al. 2013;Chanard et al. 2018); however, in many cases, the GPS time-series might contain periodic or quasi-periodic signals that cannot be explained by the loading models and are likely to be caused by phase mismodelling during the processing or due to instrumental effects. For example, Koulali & Clarke (2020) showed that GPS time-series in Antarctica might contain a quasi-annual apparent surface displacement related to snow intrusion inside the antenna through drainage holes. This signal is abruptly halted when the holes are permanently plugged, resulting in a complex pattern of a periodic signal where the annual oscillations are limited in time (Fig. 4). To demonstrate the performance of the GP method in these situations, we used a time-series of the vertical component from homogenous processing of the Antarctic network using GAMIT-GLOBK software (Herring et al. 2016). Fig. 4 shows the posterior mean of the GP fit to the uncorrected time-series data for  three GPS sites located in West Antarctica. The model fits very well the observations and reproduces the seasonal variability in the GPS position time-series. We found the maximum difference between the velocity estimated using the standard model and the GP median velocity is 0.13 mm yr −1 for BURI (Supporting Information Table  S1). Although the corrected velocity is closer to our GP median, the uncertainty is larger than the 1σ uncertainty from the GP model, which might be related to the fact that ignoring the seasonal stochastic variation results in an overestimation of the rate uncertainty. Supporting Information Table S1 summarizes the comparison between the standard model, the GP and the EMD-based approach as described in Koulali & Clarke (2020). We noted a reduction in the RMS of the GP model residuals compared to the EMD approach. We think this is a consequence of the choice of the kernel function. The GP covariance function absorbs not only the quasi-seasonal periods but also other short-term systematics. Although the EMD approach performs well in removing the snow-related spurious quasi-periodic annual signal, it relies on GRACE observations for restoring the geophysical loading signal. The advantage of the GP method here is that the covariance structure reproduces the time variability of the spurious signal without a priori removal of a surface deformation model.

Alaskan plate boundary observatory GPS
There are many GPS sites in the Plate Boundary Observatory (PBO) network that show discontinuities and quasi-periodic position signals in the time-series related to snow and ice accumulation on antennas. In some cases, the signal is quite smooth and cannot be rejected as an outlier or anomaly and it is usually retained in the time-series when estimating velocity fields. For demonstration, we selected the site AB35 which is located in Cape Yakataga, Alaska.
The time-series of this site shows high seasonal variation likely due to external rime accumulation on the antenna. Fig. 5 shows the detrended time-series of the north component and the mean model fit from the GP method. As demonstrated in the previous example the GP model reproduces very well the observed variations, whereas applying the standard model fails to fit the time-series and there are clear systematics left in the residuals (Fig. 5). The median of the marginal posterior probability distribution of the linear velocity is 15.74 +0.17 −0.18 mm yr −1 . Using the standard model we obtain a velocity of 15.33 ± 0.60. Using the MLE analysis on the residuals, we obtain a spectral index of −1.53 using the standard model and −0.38 for the GP residuals. The lack of time-variable seasonal modelling resulted in an overestimated spectral index, which led to larger uncertainties for the standard model. To validate this estimation against other methods, we used the Median Interannual Difference Adjusted for Skewness (MIDAS) method, proposed by Blewitt et al. (2016). This method was designed to provide a more robust linear velocities by mitigating the effects of seasonal signals. Using MIDAS we obtained a rate of 16.74 ± 0.25 mm yr −1 with a difference of 1 mm yr −1 with respect to the GP median, indicating that the treatment of the time variable periodic signal can lead to significant velocity differences even with long span time-series, in situations where the true signal is time varying. However, the difference between MIDAS and the GPR rate is considerably smaller than that between MIDAS and the standard model.
The GPS time-series of the site AB35 used in was derived using the GAMIT/GLOBK software. To test how our GP model performs with GPS time-series derived using a different analysis strategy, we used the time-series (AB35) from the University of Nevada, Reno (UNR) solution 3 , which is based on a precise point positioning (PPP) approach. We kept the same parametrization of the covariance matrix as described in Section 2. We found that the GP model tends to overfit the observations unless a parameter σ is added to the diagonal of the covariance matrix ('jitter' term). This means either the position uncertainties provided by the UNR solution are underestimated or the 'jitter' term absorbs the contribution from the common-mode error, which results from data processing mismodeling effects in PPP.

GRACE monthly solution
GRACE observations have been used to study a wide range of surface processes, including terrestrial water storage (TWS) and glacier dynamics. In addition to the dominant seasonal cycle, the GRACE time-series contain a significant interannual variability (e.g. Vishwakarma et al. 2021). Not accounting for this variability can lead to erroneous interpretation of the linear trend estimates of mass change. The approach proposed in this study is not restricted to GPS coordinate time-series, but can also be applied to any type of geodetic time-series. Similar to Davis et al. (2012), we applied the GP model to the equivalent water height (EWH) time-series for one of the Alaskan glaciers in the Lake Clark National Park, generated using the ANU GRACE visualization tool Darbeheshti et al. (2013). Although the noise properties and the sampling rate (monthly EWH estimates) are different from the GPS case, the GP model reproduces well the quasi-periodic signal, and the residuals do not show any obvious systematics compared to the standard model (Fig. 6a, b). Applying the time-constant model with a noise consisting of white and power-law models, we obtain a linear rate of −46.22 ± 10.04 (WRMS = 12.08), while applying the GP we infer a median rate of −41.26 +5.68 −6.22 (WRMS = 4.83). Visual inspection of the EWH time-series indicate some variability in the trend before and after 2010. However, since we are not allowing for a time-variable rate in our study, the annual amplitude may be absorbing some of the linear rate variability.

C O N C L U S I O N
For many geophysical investigations, GPS time-series are used to derive linear velocities of surface displacement. For example, Antarctic GPS vertical position time-series are fitted with trajectory models to estimate linear trends that are used for comparisons with forward GIA models (e.g. Whitehouse et al. 2012;Ivins et al. 2013), or to constrain GIA inverse solutions (e.g. Martín-Español et al. 2016) as well as inferring mantle viscosity (Barletta et al. 2018). However, isolating the so-called secular linear velocities is not a trivial task due to the presence of nonlinearities in the timeseries, including signals associated with the present-day surface mass change. This problem becomes even more complicated over short periods in the presence of large interannual variability. Therefore, it is important to consider those time-dependent signals to be able to minimize the bias.
We have proposed a method for modelling quasi-periodic signals in geodetic time-series based on GP, as an alternative to deterministic approaches. This method allows for better modelling of the systematics in GPS time-series in the presence of correlated noise. Our proposed method is not restricted to daily GPS time-series, but it applies to a different geodetic time-series with different sampling rates.
Fitting GPS time-series in areas where changes of groundwater withdrawals are significant can be a challenging task, given the time-series often show large variations in vertical surface displacement as well as highly variable annual cycle amplitudes due to substantial alternations of drought and wet periods. GP methods are a good alternative to traditional methods to mitigate the secular rate estimation bias, especially in time-series spanning short periods as well as leading to better predictions of coordinates positions.
One of the main challenges for GP modelling is the computational cost. GPs typically scale as O(N 3 ). Thus, for large problems (>10 4 data points), matrix storage and linear system solving become prohibitively expensive. This is particularly the case for Bayesian inference that requires multiple calls of the likelihood function such as MCMC methods. Olivares-Pulido et al. (2020) used MCMC for geodetic time-series trajectory modelling including noise properties. They showed for GPS time-series with 10 3 data points the computational time required is around 10 4 s CPU, which is an order of magnitude slower than using the MLE approach used in CATS software. In our study, we chose to use the CELERITE package for computational efficiency. The CPU time required for a time-series with 1278 data points is 20.03 s on an Intel i5-8350U 3.600GHz laptop. This computation gain comes at the price of the limited choice of covariance functions. The usage of more complicated covariance functions will then require using non-CELERITE functions and therefore the computational burden will increase. Another limitation of our method is related to offsets (change-point) in the time-series. In this study, we estimate the offsets as parameters of the mean function in the GP model at known epochs. However, in real case GPS time-series, offsets due to instrumental changes are not well documented by operators. A more convenient GP model for GPS time-series will include the description of multiple sudden changes in the positions in the covariance function. We leave for future investigations the usage of non-CELERITE kernel covariance functions, including other aspects of GPS time-series modelling such as post-seismic deformation and other transient deformation signals.

A C K N O W L E D G E M E N T S
Figures in this paper were produced using the Matplotlib python package and the corner (Fig. 2) was made using the corner package (https://github.com/dfm/corner.py). GPS data used in this study were obtained from UNAVCO (https://www.unavco.org) and from the International GNSS Service (IGS). AK and PJC are supported by UK Natural Environment Research Council grants NE/K004085/1 Downloaded from https://academic.oup.com/gji/article/226/3/1705/6253216 by guest on 30 July 2021 Figure 6. The equivalent water height (EWH) from the GRGS GRACE monthly solution (black) along with the mean GP model (orange) and the fit using a constant annual model. The residuals are shown for each model. and NE/R002029/1. We thank the associate editor Andrew Valentine and two anonymous reviewers for their helpful and constructive reviews.