Skip to content
BY 4.0 license Open Access Published by De Gruyter March 31, 2022

Stability analysis of the Iraqi GNSS stations

Sattar Isawi, Harald Schuh ORCID logo, Benjamin Männel ORCID logo and Pierre Sakic ORCID logo


The Iraqi GNSS network was installed in 2005 with help from the USA and UK. The network consists of seven GNSS stations distributed across Iraq. The network GNSS data have been comprehensively analyzed in this study; this, in turn, allowed us to assess the impact of various geophysical phenomena (e. g., tectonic plate motion and Earthquakes) on its positional accuracy, stability, and validity over time. We processed daily GPS data, spanning over more than five years. The Earth Parameter and Orbit System software (EPOS.P8), developed by the German Geoscience Research Center (GFZ), was used for data processing by adopting the Precise Point Positioning (PPP) strategy. The stacked time series of stations coordinates was analyzed after estimating all modeled parameters of deterministic and stochastic parts using the least-squares technique. The study confirmed a slight impact of the recent M 7.3 Earthquake on the Iraqi GNSS stations and concluded that the stations are stable over the study period (2013 up to 2018) and that the GNSS stations represent the movement of the Arabian plate.

1 Introduction

Figure 1 
The spatial distribution of the Iraqi GNSS stations. The black dashed lines are the boundary of tectonic plates. The star symbol represents the epicentre location of the M 7.3 Earthquake that hits Iraq – Iran border on November 12, 2017.

Figure 1

The spatial distribution of the Iraqi GNSS stations. The black dashed lines are the boundary of tectonic plates. The star symbol represents the epicentre location of the M 7.3 Earthquake that hits Iraq – Iran border on November 12, 2017.

The crucial part of national infrastructure is establishing an accurate spatial reference system that can be realized by a reliable geodetic network of permanent control stations. These stations are referenced to a well-defined coordinate system to be the backbone of the nation’s geospatial system (e. g., geodatabase and mapping systems).

Iraq followed this approach after the war in 2003 because reconstruction demanded implementing major engineering projects. Moreover, the implementation of these projects needs a stable and reliable geodetic basis. Therefore, in 2005, Iraq established a new geodetic reference system called the Iraqi Geospatial Reference System (IGRS). Since the establishment of IGRS, all new Iraqi engineering projects and mapping systems were based on it, specifically Oil & Gas, Transportation, Electrics, and Agricultural projects. IGRS realized as a coincidence with ITRF2000 at epoch 1997.0 but anchored to the Arabian plate by five permanent GNSS stations [32]. In principle, IGRS represents a snapshot of the Arabian plate’s dynamic on which Iraq lies. Furthermore, IGRS is a dynamic reference frame that continuously moves with the Arabian plate motion [10]. Besides, Iraq locates at the boundary between two different tectonic plates, Arabian and Eurasian plates (Fig. 1). This area is deformable and continuously unstable, as shown by several Earthquakes, for example M 7.3 Iraq – Iran border earthquake on November 12, 2017 [26]. Quantifying the influence of such geophysical processes on the stability of the Iraqi permanent GNSS stations necessitates an in-depth and comprehensive study because of its significant impact on the present and future geodetic surveying in Iraq.

The remarkable development of space geodetic techniques (e. g., Global Navigation Satellite Systems GNSS) made it easy to obtain continuous 3D positions for points on Earth’s surface for very long periods of up to many years. The obtained positions are characterized by sufficient accuracy to monitor the GNSS station’s antenna position change over time. This change indirectly refers to various physical processes in the Earth system (e. g., tectonic plate motion).

Based on this proposition, several researchers studied the movement of the Arabian plate by adopting space geodetic techniques. For example, [3] studied the Arabian plate motion by estimating the Riyadh SLR station’s position and velocity. Laser-ranging observations of 14 years from 20 SLR stations to LAGEOS-1 and LAGEOS-2 have been analyzed in this study. They concluded that the Arabian plate is countering a north-east motion with an annual linear rate of 42.9 mm/year. [11] processed ten years of daily GPS data from two Iraqi GNSS stations. ISER, located in Erbil, and ISNA, located in Najaf. The GITSA analysis software was used in this study to build up and analyze the time series of station coordinates. The study concluded that the estimated velocity vectors for the analyzed stations were 38 mm/year and 40 mm/year, respectively.

The estimation of linear rate (mainly driven by plate motion) is the main focus of these studies. However, this estimation needs an improvement to the adopted functional model. The periodic effects (e. g., annual and semi-annual deformation), sudden effects (e. g., offset due to antenna change or co-seismic deformation), or even post-seismic deformation have to be modeled and extracted from the time series signal [4]. Although all common deterministic parts would be modeled for the time series, different noise sources are still present in its signal due to satellite orbit error, atmospheric delays, and clock instability. On the other hand, the assumption that the only white noise (uncorrelated noise) contributes to the noise model will lead to underestimated uncertainty of linear velocities by factors from three to six times [33], or even five to 11 times [20].

The Iraqi GNSS stations will be extensively analyzed in this study to assess the impact of various geophysical processes (e. g., tectonic plate motion and Earthquakes) on their accuracy, stability, and validity over the study period. The determination of coordinate time series characteristics (trend, annual signals, offsets, and spectra) will be further described and discussed. Besides, this study will investigate the potential contribution of different noise types (e. g., power-law noise) to the noise model. A Matlab code comprises functions written to implement this study in both the frequency and time domain.

2 Data

Figure 2 
GPS antennas attached to the monuments. ISER, ISKU (top) and ISBA, ISNA (middle), and ZAXO station (bottom). Source NGS at

Figure 2

GPS antennas attached to the monuments. ISER, ISKU (top) and ISBA, ISNA (middle), and ZAXO station (bottom). Source NGS at

Table 1

Log files information of selected GNSS Stations.

City Najaf Kut Baghdad Erbil Dohuk
Ant. Type TRM57971 TRM57971 TRM41249 TRM57971 TPSCR.G5
Monum. Desc. Tribr. adap. Trib. adap. Trib. adap. Trib. adap. Trib. adap.

Five continuously operating GNSS reference stations have been selected out of seven Iraqi GNSS stations for this study. The selected GNSS stations are ISNA, ISKU, ISBA, ISER, and ZAXO. ISBS and ISSD stations were excluded because both of them were non-operational during the investigated period. The spatial distribution of the Iraqi GNSS stations is shown in Fig. (1). All selected GNSS stations were equipped with dual-frequency GNSS receivers as well as geodetic antenna types. However, all of them have non-geodetic monuments because of the antennas mounted on top of buildings (Fig. 2). Station ISBA is part of the International GNSS Service (IGS) network [12], whereas the remaining stations are a part of the National Geodetic Survey (NGS) global network.

Figure 3 
Number of daily observation files per station.

Figure 3

Number of daily observation files per station.

The daily observation (raw data) of all GNSS stations were downloaded from the NGS public archive at The number of daily observation files for each station is shown in Fig. (3). The GNSS station’s meta-data (log files) needed to track the station history were also obtained from the NGS database system (Table 1). The data included a comprehensive record of station firmware and hardware changes (e. g., GNSS receiver and antenna type). A pre knowledge in station history is required in this study because the change in station hardware and software could lead to offsets in the position time series [27].

2.1 GPS daily solutions

The software package EPOS.P8, developed at GFZ, was used to process the GPS data. Satellite orbits, clock corrections, and Earth Orientation Parameters (EOP) were provided by a GFZ internal reprocessing based on the operational final solution [18]. The products were given in the IGS14 reference frame. Our final solution of the GPS daily coordinates for each station was estimated following the Precise-Point-Positioning (PPP) methodology adopted in GFZ AC. The processing details (Table 11) and the flow diagram (Fig. 8) of PPP processing are shown in Appendix A.

3 Methods

Firstly, the 3-D Cartesian coordinates X i , Y i , and Z i obtained from GNSS data processing were transformed into topocentric North, East, and Up. The coordinate transformation from the right-hand system Earth Centered Earth Fixed (ECEF) into the topocentric left-hand system has been implemented according to [24]. The variance-covariance matrix SXYZ of GNSS daily solutions have also been transformed into the topocentric system by applying a variance-covariance propagation law [25].

3.1 Functional model

The time series for individual topocentric components (North, East, and Up) was built up based on the GNSS daily solutions y t i at uniformly spaced time t i ( i = 1 , ... , n, where n is the number of daily solutions). The stacked time series have been modeled (approximated) as the sum of deterministic parts, by considering the initial site position x 0 , and linear rate v x (interpreted as the GNSS station velocity), as well as the periodic signals of annual and semi-annual motion. The functional model reads:

(1) y t i = x 0 + v x t i + a sin ( 2 π t i ) + b cos ( 2 π t i ) + c sin ( 4 π t i ) + d cos ( 4 π t i ) + ε i

with a and b are the coefficients of the annual periodic motion, while c and d are the coefficients of semi annual motion. When there are in addition m offsets (jumps) in the coordinate time series, the functional model can be expanded to [6]:

(2) y t i = x 0 + v x t i + a sin ( 2 π t i ) + b cos ( 2 π t i ) + c sin ( 4 π t i ) + d cos ( 4 π t i ) + j = 1 m f j h ( t i t f j ) + ε i

where f j is the offset’s magnitude, t f j is the offset’s epoch with h denoting the Heaviside function, and ε is the model error. With apriori known offset epochs t f j , the linear model of observation equations in Eq. (2) is solved using the least-squares method:

(3) x ˆ = A T S y y 1 A 1 A T S y y 1 y ,

where the vector of adjusted unknowns x ˆ consists of the following parameters:

(4) x ˆ = x 0 ˆ v x ˆ a ˆ b ˆ c ˆ d ˆ f ˆ T ,

S y y is the variance-covariance matrix of the observation vector y,

(5) y = [ y t 1 , , y t n ] T ,

and, the design matrix A is given as:

(6) A = 1 t 1 sin ( 2 π t 1 ) cos ( 2 π t 1 ) sin ( 4 π t 1 ) cos ( 4 π t 1 ) h ( t 1 t f n ) 1 t n sin ( 2 π t n ) cos ( 2 π t n ) sin ( 4 π t n ) cos ( 4 π t n ) h ( t n t f m )

3.2 Outlier detection and removal

Outliers have been detected and removed from each coordinate time series during the preliminary model fitting. The median and interquartile range approach [21] were implemented in this study. The median and IQR statistics are calculated for the residuals on a sliding window of one year centered on each time series element to consider the gradual decrease in the data variance associated with the increase in the time series length. Outliers were defined as having:

(7) | v i m e d i a n | > 1.5 × I Q R

where I Q R defines the middle 50 % of the sample data, and v i is the i t h least square residual.

Table 2

Summary of coordinate data.

Site Data Start epoch End epoch Outliers
# # %
N 1743 2013.00 2018.99 134 7.6
E 1743 2013.00 2018.99 116 6.6
U 1743 2013.00 2018.99 145 8.3
N 1327 2013.00 2018.99 89 6.7
E 1327 2013.00 2018.99 88 6.6
U 1327 2013.00 2018.99 51 3.8
N 1591 2013.00 2018.99 70 4.3
E 1591 2013.00 2018.99 78 4.9
U 1591 2013.00 2018.99 100 6.2
N 1837 2013.00 2018.99 284 15.4
E 1837 2013.00 2018.99 267 14.5
U 1837 2013.00 2018.99 155 8.4
N 1090 2015.00 2018.99 110 10
E 1090 2015.00 2018.99 69 6.3
U 1090 2015.00 2018.99 65 5.9

Summary of removed outliers for each station time series is presented in Table (2).

4 Noise characteristics and analysis

As with many other geophysical phenomena, noise in coordinate time series can be described as a power law process of the form [2], [16]:

(8) P v ( f ) = P 0 f f 0 k

where f is the temporal frequency, P 0 and f 0 are normalizing constants, and k is the spectral index that defines a slope of the best fit line to the spectra presented in log-log space [33]. A smaller spectral index implies a more correlated process and a more relative power at lower frequencies. Therefore such a natural process is characterized by negative indices [13]. According to [17], the spectral index of different geophysical phenomena is classified into two groups. The first group is called fractional Brownian motion with spectral index 3 < k < 1. The noise processes in this interval are non-stationary including random walk with k = 2 or( P x 1 / f 2 ), while the second group termed fractional Gaussian with spectral index 1 < k < + 1, where the noise processes in this interval are stationary (independent of time) including the special case of uncorrelated white noise with k = 0 or ( P x 1 / f 0 ). The special case k = 1 or ( P x 1 / f 1 ) is known as “flicker noise”. Plenty of studies showed that the spectral index has particular importance because it is a good indicator for characterizing the noise source [14], [22]. Based on the preceding, it is necessary to study the post-fit residuals in-depth to decide which colored (correlated) noise type has to be considered in the noise model.

Two approaches are described in the literature to classify and quantify noise in the coordinate time series. The first approach’s implementation relies on time series analysis in the frequency domain (spectral analysis). In contrast, the second one relies on the time series analysis in the time domain. In this study, we restricted ourselves to the second approach, but in addition, only estimating the spectral index in the frequency domain.

4.1 Noise analysis

The time series signal and its noise do not behave in the same way when transformed and plotted as a power spectrum. This property can be of great benefit for noise analysis in the frequency domain. Based on this concept, the time series of post-fit residuals have been transformed into the frequency domain for further analysis.

Unfortunately, equipment malfunction, interruptions in the communications network, or even outlier removal lead to missing data in the coordinate time series. Therefore, a Lomb-Scargle periodogram introduced by [15] and used by [33] has been applied to solve this issue.

According to [33], the time series of post-fit residuals must be long enough and comprise data with a range of frequencies that make the power spectra well approximated. Hence, from the practical aspect, the spectral index of the power-law process can be estimated by fitting a straight line to these frequencies to avoid biasing the regression estimate with the predominance of white noise at high frequencies and decreasing the impact of outliers at lower frequencies. The functional model of power-law process in Eq. (8) can be reformulated as:

(9) P v i ( f ) = P 0 f 0 k f i k

substitution leads to

(10) W 0 = P 0 f 0 k

Taking the log of both sides yields

(11) log ( P v i ( f ) ) = log ( W 0 ) + k log ( f i )

where all parameters are defined as in Eq. (8). To estimate the slope of the fitted line (spectral index k) to the spectra at log-log space, the least-squares method was used. Therefore, Eq. (3) has been applied with the design matrix of the form:

(12) A = 1 log ( f 1 ) 1 log ( f 2 ) . . . . . . 1 log ( f N )

The spectral index has been estimated for each time series. The proper noise model was founded on a combination of white and flicker noises. Since the amplitudes of both noises have to be known in the adopted noise model, the contribution for each noise has been estimated in this study using the Maximum likelihood method presented in [30]. In this approach, the amplitudes (variances) were estimated by explicitly maximizing the likelihood function for the total variance σ 2

(13) σ 2 = v ˆ T C 1 v ˆ N

where C is the covariance matrix of post-fit residuals, N is the number of daily solutions, and v ˆ is the vector of post-fit residuals. According to [30], the covariance matrix C can be decomposed into two components with amplitudes σ w 2 and σ k 2 being the variance of white and power-law noises respectively such that:

(14) C = σ w 2 I w + σ k 2 J k

where Iw is the identity matrix, and Jk is the covariance matrix of power-law noise, and both of them have the size of N × N. The elements of matrix Jk with k = 1 (flicker noise) have been approximated in this study according to [33]:

(15) J -1 3 4 2 × 2 t i = t j 3 4 2 × 2 log ( t i t j ) log ( 2 ) + 2 12 t i t j

The approximated covariance matrix in Eq. (15) is not exactly the same as that derived by the method of fractional difference, introduced and used by [9], [28], and [30]. The only difference is in the scaling [29].

5 Results and discussion

High precision coordinate time series have been produced after processing GNSS data and cleaning for outliers. The most common deterministic and stochastic parts were modeled and described in Eq. (2). Various model parameters were estimated and assessed. In this section, we provide further details and discuss all estimated parameters.

5.1 Stochastic part

It is necessary to investigate the post-fit residuals further to build up a proper noise model. Hence, real uncertainties can be assigned to the estimated long-term rate (velocity) and other parameters in the underlying functional model. A typical example of estimated power spectral density and its associated spectral index is given for the station ISKU. The station observations span over more than five years of GNSS data as presented in Fig. (4). At this figure, the noise power is almost constant for up to cross-over frequency at (1/12  day) which implies the predominant white noise at high frequencies. In contrast, the power increases from a cross-over frequency to the left, indicating the predominant time-dependent noise at low frequencies.

Figure 4 
The Lomb-Scargle periodogram for the station ISKU. The solid lines (magenta) are the best fit in log-log space from frequencies less than the cross-over frequency 1/(14 days).

Figure 4

The Lomb-Scargle periodogram for the station ISKU. The solid lines (magenta) are the best fit in log-log space from frequencies less than the cross-over frequency 1/(14 days).

Figure 5 
Estimated spectral indices.

Figure 5

Estimated spectral indices.

Furthermore, due to errors in the solid Earth tide model, or in the indirect tidal errors in the satellite orbits [1], a noise peak at the frequency (1/13.6  day) is observed in this site signal, and at most analyzed time series as well. However, the estimation of spectral indices was done for frequencies more minor than (1/14  day) to prevent biasing the estimation with this peak and exclude the predominant white noise.

Table 3

Estimated white + flicker noise amplitudes.

Site Spectral Index W.Noise ± σ (mm) % F.Noise ± σ (mm) %
N −0.46 ± 0.06 0.41 ± 0.03 5 3.16 ± 0.23 95
E −0.79 ± 0.06 0.65 ± 0.16 6 4.32 ± 0.02 94
U −0.68 ± 0.05 1.45 ± 0.02 8 8.45 ± 0.08 92
N −0.13 ± 0.06 0.19 ± 0.02 1 3.31 ± 0.29 99
E −0.46 ± 0.05 0.31 ± 0.24 2 4.09 ± 0.02 98
U −0.47 ± 0.06 1.14 ± 0.01 5 8.69 ± 0.10 95
N −0.41 ± 0.06 0.58 ± 0.03 6 3.85 ± 0.19 94
E −0.59 ± 0.06 0.81 ± 0.15 8 4.85 ± 0.03 92
U −0.54 ± 0.06 1.29 ± 0.01 6 9.21 ± 0.08 94
N −0.51 ± 0.05 0.76 ± 0.03 10 4.01 ± 0.17 90
E −0.73 ± 0.05 0.69 ± 0.15 6 4.63 ± 0.02 94
U −0.54 ± 0.05 1.03 ± 0.01 4 9.10 ± 0.08 96
N −0.24 ± 0.05 0.25 ± 0.02 1 3.88 ± 0.32 99
E −0.51 ± 0.06 0.58 ± 0.25 5 4.65 ± 0.03 95
U −0.33 ± 0.06 1.36 ± 0.02 6 9.22 ± 0.12 94

The spectral analysis shows that the estimated spectral indices (Fig. 5) range from −0.20±0.05 to −0.50±0.04 for the north direction and −0.48±0.04 to −0.59±0.04 for the east direction, while the indices for up direction range from −0.37±0.08 to −0.57±0.03. Furthermore, the estimated spectral indices for all time series are ranging within fractional Gaussian part ( 1 < k < 0). This means the power-law noise being the main contributor to the character of stochastic part, confirming the findings of [33], [20], and [30]. Hence, the noise model of each time series was determined accordingly.

Further to preceding spectral analysis, the adopted noise model identified previously was analyzed. This analysis included quantifying the contribution of white and flicker noises by estimating the amplitudes (variances) in the time domain. The analysis shows a predominant contribution for flicker noise to the total noise model. This contribution ranges from 90 % to 99 %, while the white noise contribution ranges from 1 % to 10 %. The initial findings suggested to perform the data fitting once more for all coordinate time series, but in this case, with a presumed stochastic model mainly based on the estimated contribution for each noise type. The estimated amplitudes are provided in Table (3). This Table shows that the estimated amplitudes of flicker noise are very consistent for all coordinate time series. The same also applies to white noise, indicating that a similar noise type contaminates all analyzed coordinate time series.

Figure 6 
Estimated offsets due to M 7.3 Earthquake.

Figure 6

Estimated offsets due to M 7.3 Earthquake.

5.2 Regression analysis

The completion of the noise analysis created a solid ground to implement the regression analysis. This, in turn, allowed us to adopt a realistic stochastic model to proceed with the weighted least-squares estimation for the data model in Eq. (2). Various model parameters were estimated and assessed in this analysis, for example, the time series jump (mainly driven by Earthquakes), trend (plate motion), and the annual and semi-annual variation of coordinate time series. In this section, we discuss in detail the outcomes of the regression analysis we have performed.

5.2.1 Discontinuities due to Earthquake

We evaluated the impact of the recent M 7.3 Earthquake on the stability of the network stations. The evaluation has been done simultaneously with testing the slight potential effect of receiver hardware and firmware change on the coordinate time series continuity. The findings of this evaluation are presented in Fig. (6). This figure shows that all stations displaced in the up direction except the station ISER which is only displaced in the north direction. ZAXO station was relatively more influenced (4.5 mm offset at up direction) than the remaining stations whose offsets ranged from 2.4  mm to 3.40  mm at up direction. Details of all observed offsets caused by Earthquake and firmware change are provided in Table (4).

Table 4

Observed Offsets.

Site Axis Offset (mm) Epoch Caused by
ISBA E 1.70 ± 0.25 2017.8657 M 7.3 earthquake
ISBA U 2.50 ± 0.47 2017.8657 M 7.3 earthquake
ISBA U −3.50 ± 0.65 2015.2054 firmware change
ISNA E −1.01 ± 0.43 2017.7534 firmware change
ISNA U 2.40 ± 0.95 2017.8657 M 7.3 earthquake
ISKU N −1.20 ± 0.30 2017.8657 M 7.3 earthquake
ISKU U 3.40 ± 0.71 2017.8657 M 7.3 earthquake
ISKU U −1.23 ± 0.60 2015.1643 firmware change
ISKU U 1.97 ± 0.68 2016.9890 firmware change
ISER N −1.10 ± 0.25 2017.8657 M 7.3 earthquake
ISER N 2.40 ± 0.23 2014.0958 unknown
ISER N 1.30 ± 0.21 2014.6876 firmware change
ISER E 1.40 ± 0.27 2014.0958 unknown
ISER E −1.20 ± 0.55 2014.6876 firmware change
ZAXO E 1.30 ± 0.35 2016.1041 firmware change
ZAXO U −1.70 ± 0.69 2016.1041 firmware change
ZAXO U 4.50 ± 0.79 2017.8657 M 7.3 earthquake

Table 5

Estimated staions velocities.

Site Velo. (mm/yr) σ (mm/yr) Velo. (mm/yr) σ (mm/yr)
White+Flicker White-only
N 27.67 ±0.13 27.79 ±0.03
E 24.11 ±0.16 24.07 ±0.04
U −0.92 ±0.24 −1.25 ±0.08
N 27.27 ±0.10 27.20 ±0.04
E 24.85 ±0.12 24.73 ±0.05
U −1.98 ±0.23 −1.93 ±0.09
N 26.18 ±0.15 26.25 ±0.04
E 23.26 ±0.18 23.29 ±0.05
U −3.30 ±0.23 −3.51 ±0.09
N 25.03 ±0.18 24.96 ±0.04
E 22.14 ±0.17 21.87 ±0.05
U −4.92 ±0.20 −5.05 ±0.08
N 26.27 ±0.14 26.82 ±0.08
E 22.82 ±0.21 23.00 ±0.09
U −1.39 ±0.32 −2.02 ±0.18

5.2.2 Station velocities

The outcomes of the regression analysis show very similar horizontal velocities for all stations (Table 5). It is an expected result as all stations are close to each other within the rigid part of the Arabian plate. The estimated velocities after removing all observed offsets (jumps) range from 25.0 ± 0.2 m m / y e a r to 27.7 ± 0.1 m m / y e a r for the north direction and 22.1 ± 0.2 m m / y e a r to 24.9 ± 0.1 m m / y e a r for the east direction. The velocity vectors depicted in Fig. (7) shows that all stations experienced a north-east movement, which mainly represents a snapshot of the Arabian plate motion in Iraq. In contrast, our study found that the vertical velocities 4.9 ± 0.2 m m / y e a r and 3.3 ± 0.2 m m / y e a r for ISER and ISKU respectively (see Appendix B and C), are remarkable vertical velocities and worthy of attention and assessment. The most realistic explanation to the observed velocities might be the effect driven by non-tectonic signals. This affection can cause deformation, for example, groundwater withdrawal produced long-term subsidence [5]. The less likely but interesting explanation is that all network stations have a non-geodetic monument because the antennas were mounted on buildings (Fig. 2). Consequently, the buildings themselves could be subject to deformation. On the other hand, we estimated the velocity rates for all stations two times by applying two different stochastic models to emphasize the importance of the implemented stochastic model in this study. First, the velocity rate was estimated using only a white noise-based stochastic model, while at the second time, a stochastic model comprised white noise plus power-law noise (flicker noise). Table (5) shows the estimated velocities and their uncertainties in both cases. As expected, the resulting uncertainties for case one (white noise-based stochastic model) are much smaller than the second case. The underestimation by almost one order of magnitude confirms the findings of [33], and [21].

Figure 7 
Locational map of the M 7.3 Earthquake epicenter (gray star) and the depicted velocity vectors (black arrows) of the Iraqi GNSS stations.

Figure 7

Locational map of the M 7.3 Earthquake epicenter (gray star) and the depicted velocity vectors (black arrows) of the Iraqi GNSS stations.

5.2.3 Annual variations

Annual and semi-annual variations have been observed in all coordinate time series, at both horizontal and vertical components. The annual variations range from 0.5 ± 0.1 m m to 2.5 ± 0.1 m m for the north direction and 0.2 ± 0.1 m m to 3.9 ± 0.1 m m for the east direction, while in the up direction the situation is quite different. The annual variations at vertical components range from 7.8 ± 0.2 m m to 12.8 ± 0.1 m m. Although these results are in good agreement with each other, the loading corrections, as suggested for example by [19] were not applied in processing GNSS data. Thus, the height variations are much larger than the variations in the horizontal components. Furthermore, the local hydrological loading causes height variations of (10–30  mm) with a predominant annual period [8]. For example, the Tigris and Euphrates rivers and the several dams distributed at middle and north of Iraq. These rivers and dams are usually full in the winter and are relatively empty in the summer season. Besides, the thermal expansion of the station monuments could contribute to this oscillation [23]. The stations monuments are exposed to a significant temperature variation due to the temperature change in Iraqi’s summer and winter season. The average daily temperature in winter reaches 6 °C, while 48 °C in summer [31]. Hence, the station monument may experience a thermal expansion because the GNSS antennas are mounted on galvanized iron rods at top of the buildings (Fig. 2). Details of estimated annual and semi-annual variations are presented in Tables (6) and (7), respectively. The post-fit residuals obtained after removing the annual variations from all coordinate time series are presented in Table (8). This table shows very similar residuals for all analyzed time series, indicating a similar noise type.

Table 6

Estimated annual signals.

Site Amplitude (mm) Uncertainty (mm) Phase (deg) Uncertainty (deg)
N 1.49 ±0.07 −118.32 ±0.002
E 2.02 ±0.09 −080.95 ±0.002
U 12.85 ±0.05 −112.74 ±0.006
N 0.88 ±0.08 −163.57 ±0.006
E 0.24 ±0.10 021.99 ±0.007
U 10.98 ±0.20 −112.11 ±0.008
N 2.46 ±0.08 082.00 ±0.001
E 1.00 ±0.11 −089.51 ±0.001
U 10.14 ±0.20 −117.66 ±0.008
N 0.50 ±0.08 −141.63 ±0.006
E 2.61 ±0.10 −085.92 ±0.001
U 7.72 ±0.18 −108.17 ±0.006
N 0.73 ±0.10 −069.68 ±0.004
E 3.86 ±0.11 104.46 ±0.003
U 12.39 ±0.24 −134.17 ±0.015

Table 7

Estimated semi-annual signals.

Site Amplitude (mm) Uncertainty (mm) Phase (deg) Uncertainty (deg)
N 1.04 ±0.05 141.76 ±0.002
E 1.19 ±0.05 130.54 ±0.001
U 0.95 ±0.01 117.74 ±0.003
N 0.81 ±0.07 130.18 ±0.015
E 0.81 ±0.33 107.73 ±0.006
U 0.74 ±0.01 047.50 ±0.005
N 1.17 ±0.04 172.98 ±0.001
E 0.74 ±0.08 155.97 ±0.001
U 1.40 ±0.03 061.36 ±0.004
N 0.58 ±0.09 178.21 ±0.013
E 0.80 ±0.03 151.74 ±0.001
U 1.98 ±0.04 085.68 ±0.001
N 0.75 ±0.11 113.82 ±0.001
E 0.80 ±0.02 155.34 ±0.003
U 2.40 ±0.04 048.98 ±0.014

Table 8

Post-fit residuals scatter.

Site N (mm) E (mm) U (mm)
ISBA 1.90 2.60 5.10
ISNA 1.90 2.40 5.10
ISKU 2.30 2.90 5.40
ISER 2.40 2.70 5.30
ZAXO 2.20 2.70 5.50

Table 9

Summary of first comparsion (GNSS vs. SLR).

Study Site Estimated Velocity Velocity. Vec. (mm/yr) Azimuth Technique

Riyadh 29.10 31.60 1.90 42.90 N47°24´E SLR
⋆⋆ ISBA 27.67 24.11 −0.92 36.70 N41°03´E GNSS
ISNA 27.27 24.85 −1.98 36.70 N42°20´E GNSS
ISKU 26.18 23.26 −3.30 35.02 N41°36´E GNSS
ISER 25.03 22.14 −4.92 33.41 N41°29´E GNSS
ZAXO 26.27 22.82 −1.39 34.79 N40°58´E GNSS

  1. ⋆ Alothman and Schillak [3] ⋆⋆ This study

Table 10

Summary of second comparsion.

Study Site Velocity Vector (mm/yr) Azimuth Space Technique
ISER 38.00 N40°48´E GNSS
ISNA 40.00 N45°06´E GNSS
⋆⋆ ISER 33.41 N41°29´E GNSS
ISNA 36.70 N42°20´E GNSS

  1. ⋆ Jasim et al. [11] ⋆⋆ This study

6 Validating the results

The estimated velocities of the GNSS stations in this study were validated with two independent studies. However, the two studies did not account for the discontinuities (offsets) and the annual variation of the coordinate time series. The velocity vectors estimated by [3] using Satellite Laser Ranging for the SLR-station at Kingdom of Saudi Arabia, also located at the stable part of the Arabian plate, agrees very well with the velocity rates observed in this study (Table 9). The second validation (Table 10) shows, that the velocities of the current study are 4 mm/year less than those estimated by [11] that was obtained from processing ten years of daily observation data for two GNSS stations in mid and north of Iraq.

7 Conclusion

We assessed the stability of Iraq’s core GNSS stations based on the processing and interpretation of the GNSS daily data from 2013 up to 2018, derived for five permanent GNSS stations across Iraq. Determining coordinate time series characteristics (trend, annual signals, discontinuities, and spectra) were comprehensively discussed.

The spectral analysis showed that the estimated spectral indices for all time series fall within the fractional Gaussian part ( 1 < k < 0). This, in turn, refers to the power-law noise being the main contributor to the stochastic part. Besides, the estimated annual signals with amplitudes of around 10 mm were well observed in the up component for all time series. The signals agreed with each other and reflected the impact of different geophysical processes in the Earth system. The study concluded, that the source of such annual signals mainly refers to the near-surface hydrologic loading and the seasonal temperature variation.

The study shows a slight impact of the recent M 7.3 Earthquake on the Iraqi GNSS stations. The maximum offset does not exceed five milimeters. Moreover, our analysis indicates that several offsets caused by equipment firmware change occurred; and were removed in the coordinate time series.

The determination of station velocities included a comparison between different noise assumptions. The study reveals the benefit of choosing a suited noise model (e. g., the combination of white and power-law noises) on the derived velocities. The estimated stations velocities range from 25 ± 0.1 m m / y e a r to 27.6 ± 0.1 m m / y e a r for north direction and 22.1 ± 0.1 m m / y e a r to 24.8 ± 0.1 m m / y e a r for east direction. Subsidence rates of −4.9  mm/year and −3.3  mm/year are observed for stations ISER and ISKU, respectively.

A comparison was made with two former studies to validate the results of this study. The validation showed very good agreement with the results of [3], [11], and the study confirmed that the Iraqi GNSS stations are stable over the study period and that the stations represent the Arabian plate’s motion.

Appendix A PPP processing informations

Table 11

Basic models and parameters of PPP processing.

Modeling and a priori information
Observations ionosphere free linear combination formed by undifferenced GPS observation, sampled with 300  s, elevation cut-off 5°
A priori products GFZ SP3 orbits, clock corrections, Earth orientation parameters [18]
Tropospheric correction empirical model GPT2 mapped with GMF [7]
Ionospheric correction first-order effect considered with ionosphere-free linear combination
GNSS phase center igs14.atx
Estimated parameters
Coordinates estimated daily in IGS14 reference frame introduced via the orbit products
Troposphere 25 zenith delays per day; two gradient pairs per station and day
Receiver clock pre-eleminated every epoch
GNSS ambiguities estimated, without ambiguity fixing
Figure 8 
Flow diagram of PPP processing via EPOS.P8 software.

Figure 8

Flow diagram of PPP processing via EPOS.P8 software.

Appendix B Position time series of site ISKU

Figure 9 
The Top, middle, and bottom are the position time series for north, east, and up respectively. The daily solutions (black dots) of the site ISKU are modeled (red lines) with a linear trend. The blue lines refer to the removed related offsets.

Figure 9

The Top, middle, and bottom are the position time series for north, east, and up respectively. The daily solutions (black dots) of the site ISKU are modeled (red lines) with a linear trend. The blue lines refer to the removed related offsets.

Appendix C Position time series of site ISER

Figure 10 
The Top, middle, and bottom are the position time series for north, east, and up respectively. The daily solutions (black dots) of the site ISER are modeled (red lines) with a linear trend. The blue lines refer to the removed related offsets.

Figure 10

The Top, middle, and bottom are the position time series for north, east, and up respectively. The daily solutions (black dots) of the site ISER are modeled (red lines) with a linear trend. The blue lines refer to the removed related offsets.


[1] Abraha, K.E., F.N. Teferle, A. Hunegnaw, and R. Dach. 2017. GNSS related periodic signals in coordinate time-series from precise point positioning. Geophysical Journal International 208(3): 1449–1464. DOI:10.1093/gji/ggw467.Search in Google Scholar

[2] Agnew, D.C. 1992. The time-domain behavior of power-law noises. Geophysical research letters 19(4): 333–336.10.1029/91GL02832Search in Google Scholar

[3] Alothman, A. and S. Schillak. 2014. Recent results for the arabian plate motion using satellite laser ranging observations of riyadh SLR station to LAGEOS-1 and LAGEOS-2 satellites. Arabian Journal For Science and Engineering 39(1): 217–226. DOI:10.1007/s13369-013-0823-7.Search in Google Scholar

[4] Altamimi, Z., P. Rebischung, L. Métivier, and X. Collilieux. 2016. ITRF2014: A new release of the International Terrestrial Reference Frame modeling nonlinear station motions. Journal of Geophysical Research: Solid Earth 121(8): 6109–6131. DOI:10.1002/2016JB013098.Search in Google Scholar

[5] Bawden, G.W., W. Thatcher, R.S. Stein, K.W. Hudnut, and G. Peltzer. 2001. Tectonic contraction across los angeles after removal of groundwater pumping effects. Nature 412(6849): 812–815. DOI:10.1038/35090558.Search in Google Scholar PubMed

[6] Bevis, M. and A. Brown. 2014. Trajectory models and reference frames for crustal motion geodesy. Journal of Geodesy 88(3): 283–311. DOI:10.1007/s00190-013-0685-5.Search in Google Scholar

[7] Böhm, J., B. Werl, and H. Schuh. 2006. Troposphere mapping functions for gps and vlbi from ecmwf operational analysis data. Journal of geophysical research 111(B2): B02406.10.1029/2005JB003629Search in Google Scholar

[8] Dill, R. and H. Dobslaw. 2013. Numerical simulations of global-scale high-resolution hydrological crustal deformations. Journal of Geophysical Research: Solid Earth 118(9): 5008–5017. DOI:10.1002/jgrb.50353.Search in Google Scholar

[9] Hosking, J. 1981. Fractional differencing. biometrika 68: 165–176. Mathematical Reviews (MathSciNet): MR614953 Zentralblatt MATH 464.10.1093/biomet/68.1.165Search in Google Scholar

[10] International Association of Oil & Gas Producers. 2015, may. Geodetic referencing in Iraq. Technical report.Search in Google Scholar

[11] Jasim, Z.N., O.Y. Alhamadani, and M.U. Mohammed. 2018. Investigation the plate motion using operating refere. International Journal of Civil Engineering 9(13).Search in Google Scholar

[12] Johnston, G., A. Riddell, and G. Hausler. 2017. The international GNSS service, In Springer handbook of global navigation satellite systems, 967–982. Springer.10.1007/978-3-319-42928-1_33Search in Google Scholar

[13] Kenyeres, A. and C. Bruyninx. 2009. Noise and periodic terms in the EPN time series, In Geodetic Reference Frames, 143–148. Springer. DOI:10.1007/978-3-642-00860-3_22.Search in Google Scholar

[14] Langbein, J., F. Wyatt, H. Johnson, D. Haman, and P. Zimmer. 1995. Improved stability of a deeply anchored geodetic monument for deformation monitoring. Geophysical Research Letters 22(24): 3533–3536. DOI:10.1029/95GL03325.Search in Google Scholar

[15] Lomb, N. and J. Scargle. 1976. Further development and properties of the spectral analysis by least-squares. Astrophysics and Space Science 39: 447.10.1007/BF00648343Search in Google Scholar

[16] Mandelbrot, B.B. 1983. The fractal geometry of nature/revised and enlarged edition. New York, WH Freeman and Co., 1983, 495 p.Search in Google Scholar

[17] Mandelbrot, B.B. and J.W. Van Ness. 1968. Fractional brownian motions, fractional noises and applications. SIAM review 10(4): 422–437. DOI:10.1137/1010093.Search in Google Scholar

[18] Männel, B., A. Brandt, T. Nischan, A. Brack, P. Sakic, and M. Bradke. 2020. GFZ final product series for the International GNSS Service (IGS). DOI:10.5880/GFZ.1.1.2020.002.Search in Google Scholar

[19] Männel, B., H. Dobslaw, R. Dill, S. Glaser, K. Balidakis, M. Thomas, and H. Schuh. 2019. Correcting surface loading at the observation level: impact on global gnss and vlbi station networks. Journal of Geodesy 93(10): 2003–2017. DOI:10.1007/s00190-019-01298-y.Search in Google Scholar

[20] Mao, A., C.G. Harrison, and T.H. Dixon. 1999. Noise in GPS coordinate time series. Journal of Geophysical Research: Solid Earth 104(B2): 2797–2816. DOI:10.1029/1998JB900033.Search in Google Scholar

[21] Nikolaidis, R. 2002, January. Observation of geodetic and seismic deformation with the Global Positioning System. Ph. D. thesis, University of California, San Diego.Search in Google Scholar

[22] Nistor, S. and A. Buda. 2014. Time series analysis in the presence of coloured noise. Journal of Applied Engineering Science 4(17): 2.Search in Google Scholar

[23] Romagnoli, C., S. Zerbini, L. Lago, B. Richter, D. Simon, F. Domenichini, C. Elmi, and M. Ghirotti. 2003. Influence of soil consolidation and thermal expansion effects on height and gravity variations. Journal of Geodynamics 35(4–5): 521–539. DOI:10.1016/S0264-3707(03)00012-7.Search in Google Scholar

[24] Seeber, G. 2003. Satellite geodesy, 2nd completely revised and extended edition. Walter de Gruyter GmbH & Co. KG 10785: 303–304.10.1515/9783110200089Search in Google Scholar

[25] Soler, T. and M. Chin 1985. On transformation of covariance matrices between local cartesian coordinate systems and commutative diagrams, In ASP-ACSM Convention, 393–406.Search in Google Scholar

[26] U.S. Geological Survey. 2017. Earthquake lists, Maps, and Statistics, accessed december 18, 2017 at.Search in Google Scholar

[27] Wanninger, L. 2009. Correction of apparent position shifts caused by GNSS antenna changes. GPS solutions 13(2): 133–139. DOI:10.1007/s10291-008-0106-z.Search in Google Scholar

[28] Williams, S. 2003. The effect of coloured noise on the uncertainties of rates estimated from geodetic time series. Journal of Geodesy 76(9–10): 483–494. DOI:10.1007/s00190-002-0283-4.Search in Google Scholar

[29] Williams, S.D. 2008. CATS: GPS coordinate time series analysis software. GPS solutions 12(2): 147–153. DOI:10.1007/s10291-007-0086-4.Search in Google Scholar

[30] Williams, S.D., Y. Bock, P. Fang, P. Jamason, R.M. Nikolaidis, L. Prawirodirdjo, M. Miller, and D.J. Johnson. 2004. Error analysis of continuous GPS position time series. Journal of Geophysical Research: Solid Earth 109(B3). DOI:10.1029/2003JB002741.Search in Google Scholar

[31] World 2021. The climate in iraq.Search in Google Scholar

[32] Yenter, M., T. Wheeler, J. Mejia, K. Joyce, and M. Cline. 2005, 04. Developing IGRS. Engineer 36: 17–18.Search in Google Scholar

[33] Zhang, J., Y. Bock, H. Johnson, P. Fang, S. Williams, J. Genrich, S. Wdowinski, and J. Behr. 1997. Southern california permanent GPS geodetic array: Error analysis of daily position estimates and site velocities. Journal of geophysical research: solid earth 102(B8): 18035–18055. DOI:10.1029/97JB01380.Search in Google Scholar

Received: 2022-01-09
Accepted: 2022-02-20
Published Online: 2022-03-31
Published in Print: 2022-07-26

© 2022 Isawi et al., published by De Gruyter

This work is licensed under the Creative Commons Attribution 4.0 International License.