Topo-Iberia CGPS network: a new 3D crustal velocity field in the Iberian Peninsula and Morocco based on 11 years (2008–2019)

The Topo-Iberia network covering the Spanish part of the Iberian Peninsula and Morocco was established in 2008. Following the first publication of the horizontal velocity field based on an analysis of the first 4 years (2008–2012) in 2015, a new 3D velocity field based on 11 years (2008–2019) is presented. The double-differencing technique is applied with GAMIT software, and the precise point positioning (PPP) technique is applied with GipsyX software for GPS processing to obtain loosely constrained solutions. Then, a kalman filter is used to merge both solutions using GLOBK software and produce a velocity solution. The combined velocity solution is presented as a new horizontal velocity field and the first vertical velocity field of the Topo-Iberia CGPS network. A global view of the data quality from the Topo-Iberia stations is also presented for the first time, as well as a program specifically created to convert GipsyX gdcov files to SINEX files, gdcov2sinex.sh. The results in the IGb14 reference frame are in good agreement with previous publications, except in the Betic Cordillera, where the new velocities with respect to Eurasia are generally lower and present a higher southwestward component rather than a purely westward motion. The main discrepancy is observed at the ALJI station, where the maximum velocity is reached (3 ± 0.1 mm/yr), showing west-northwestward motion. The smallest and highest residual displacements with respect to Eurasia occur in northern Iberia and Morocco, respectively. The results support the conclusions of other studies that a clockwise block rotation occurs at present in the northern branch of the Gibraltar Arc. The velocities indicate a northwest-southeast convergence rate of 3.5 ± 0.1 mm/yr with respect to Nubia. A general low subsidence (mean value of 0.5 ± 0.1 mm/yr) is observed in the Iberian Peninsula except for some isolated stations presenting uplift.


Introduction
The Topo-Iberia network was part of a research project entitled Geosciences in Iberia: Integrated studies of topography and 4-D evolution. 'Topo-Iberia' supported by the Spanish Government (Gallart 2006). Using a Global Navigation Satellite System (GNSS) technique, specifically the Global Positioning System (GPS), was one of the methodologies used to take crustal deformation measurements. As part of the project, a Topo-Iberia network composed of 26 continuously operating GPS stations (CGPS), covering the Spanish part of the Iberian Peninsula (22) and Morocco (4), was established ( Fig. 1). The Spanish Geological Survey (IGME) implemented the Topo-Iberia Information System (SiTopo) to store the data. Garate et al. (2015) presented a velocity field based on the analysis of the first 4 years of continuous GPS observations from the Topo-Iberia network. Data analyses were performed by three different groups using different scientific software packages and approaches. The final solution from each group was transformed to a common reference frame, and finally, an average velocity field from the three solutions was calculated.
At the time of this study, there are more than 11 years of data from some Topo-Iberia stations whose analysis can allow us to estimate the velocity field with greater precision. The abovementioned study on the Topo-Iberia network estimated only a horizontal (2D) velocity field, and 1 3 155 Page 2 of 15 few studies, such as Serpelloni et al. (2013), Nguyen et al. (2016), Devoti et al. (2017) and Khazaradze et al. (2019), have estimated the vertical deformation in the Iberian Peninsula. Therefore, the main aims of this study are (1) to process all available data from the Topo-Iberia CGPS network; (2) to offer a new 3D velocity field showing present-day deformation in the Iberia plate and along its boundary zone with the Nubia plate in the western Mediterranean region; and (3) to use software with different GNSS processing techniques. Some reasons for applying different software are to check for possible systematic effects due to software or to remove mistakes introduced during processing, which are difficult to detect without a second analysis for comparison. In this way, among software programs using double-differencing and precise point positioning (PPP) techniques, GAMIT version 10.71 (Herring et al. 2018) and the Jet Propulsion Laboratory's (JPL) GipsyX software version 1.7 (Bertiger et al. 2020) have been chosen, respectively, and GLOBK is used to analyze the time series and fix the reference system. For geodetic applications and static positioning, with 24-h sessions, the expected precision using GAMIT is at the level of 1-2 mm and 3-5 mm for horizontal coordinates and heights, respectively. On the other hand, the expected precision using GipsyX, PPP with ambiguity resolution, is 2 mm and 6.5 mm, also, respectively.
Although some EUREF Permanent Network (EPN) stations have been included in the new analysis, this work is not on the densification of the EPN, nor is it intended to be, nor is it related to any processing strategy that EUREF has established for its data. As indicated above, one of the objectives is to re-process Topo-Iberia GPS data with two different software, and for this purpose, the combination and  Mancilla et al. (2015) and personal communication with Sánchez-Gómez M (2022). The gray vector represents the convergence direction and rate between the plates based on the GSRM v2.1 model. The shaded-relief image of the Earth illumination from the south is at 1 arc minute resampled from the SRTM15 + V topography database (Tozer et al. 2019) plotted with Generic Mapping Tools (GMT) software version 5.4.5 (Wessel et al. 2013) processing strategy established by Herring et al. (2016) is followed and the International Terrestrial Reference Frame (IGb14) is used. Considering the complexity of the process, details and explanations of some strategy steps are provided to help GAMIT/GLOBK and GipsyX users. Due to space limitations, much of the detail is provided in supplementary material.

Tectonic setting
In terms of tectonic regime, the Iberian Peninsula can be divided into three main domains: (1) the Pyrenean-Cantabrian system (in the north) (Fig. 1) is a collisional orogen with double vergence, resulting from the convergence of the Iberia and Europe plates between the late Cretaceous and the early Miocene (Choukroune 1992); (2) the Central region includes deformed regions, mainly of Neogene age, in the plate interior, such as the Central System and the Iberian Chain, together with their associated sedimentary basins; and (3) the Betic-Rif system (in the south) includes the northwest-southeast Eurasia-Nubia convergent plate boundary in the westernmost Mediterranean originating in the Cretaceous (Dewey et al. 1989). In a Eurasian-fixed reference frame, the largest horizontal deformation is located in this third domain. The plates converge in this region at a rate of 4.9 mm/yr (represented in the map figures with a gray vector) according to the GSRM v2.1 plate motion model (Kreemer et al. 2014), calculated using the Plate Motion Calculator of the Geodetic Facility for the Advancement of Geoscience (GAGE) operated by the EarthScope Consortium (https:// www. unavco. org/ softw are/ geode tic-utili ties/ plate-motion-calcu lator/ plate-motion-calcu lator. html).
Numerous studies have been carried out on the deformation along the Eurasia-Nubia plate boundary zone in this region, such as Koulali et al. (2011), who established westsouthwest movements with respect to Eurasia (∼2-3 mm/ yr) in the southwestern part of the Betic Cordillera. The Betic-Rif deformation zone is controlled by slab rupture and Alboran Domain collision (Mancilla et al. 2015). In the eastern Betic Cordillera, within the complex general context of extension and westward displacements, there are moderate active compressive tectonics (Palano et al. 2013;Galindo-Zaldivar et al. 2015), where part of the convergence between Nubia and Eurasia is accommodated by shortening (Borque et al. 2019). In the western Betic Cordillera, the development of the Gibraltar Arc, the tectonic belt connecting the Betic and Rif Cordilleras surrounding the Alboran Sea, continues at present (Gonzalez-Castillo et al. 2015), and the displacement pattern is in agreement with a clockwise rotation of the tectonic units in the northern branch of the Gibraltar Arc. On the other hand, Galindo-Zaldivar et al. (2022) established the present-day interaction of the rollback tectonics of the Betic-Gibraltar-Rif slab in the western part of the Gibraltar Arc, with indentation tectonics affecting the eastern and southern areas. The Betic-Rif deformation zone and hypothetical Eurasia-Nubia plate boundary, both based on Mancilla et al. (2015) and personal communication with Sánchez-Gómez M (2022), are denoted in the map figures with a magenta and black line, respectively.

GPS data
At the time of this study, all available data from the Topo-Iberia CGPS network cover the interval 2008-2019. Over time, some stations have been dismantled or have suffered incidents that remain unresolved, and as a consequence, the continuity of data spans varies significantly between stations ( Fig. 2 top). Blewitt and Lavallée (2002) established the relationship between the data span and the effect of annual signals in velocity estimation. The recommended standard minimum data span for velocity solutions intended for tectonic interpretation or for establishing a reference frame is 2.5 years, and the velocity bias is a minimum at integer-plus-half years. Although the available data span is long, these recommendations have been taken into account. Therefore, the available span is shortened as necessary, with a priority of eliminating the periods of malfunction ( Fig. 3) with the help of a previous data quality check carried out using TEQC software (Estey and Meertens 1999).  complete observations recorded at each Topo-Iberia station is in the range of 21,640-24,489. The number of stations as a function of the selected spans is shown in Fig. 2 (bottom), highlighting the largest number of sites over 9.5 years of observation. All stations together form a dataset spanning from 2008-04-09 to 2019-10-29, resulting in a total of 11.6 years. Details on available and processed data spans are shown in Table S1 of Supplement 1 (SI_1).
Continuously operating reference stations (CORS) from the EPN are also included in this analysis to tie the velocity solution into a global reference frame. A set of ten or more well-distributed sites with small uncertainties is sufficient to obtain a robust frame realization (Herring et al. 2018). To select the best reference stations, an interactive tool provided by EUREF (http:// www. epncb. oma. be/_ produ ctsse rvices/ Refer enceF rame/) was used. Thirteen stations were selected between defined stations as recommended to be reference stations for the given span: ALAC, ALME, BELL, CACE, CANT, CASC, HERS, MELI, RABT, SALA, VALE, VIGO, and YEBE. All of them, except BELL, CACE, CASC, and MELI, are classified by EUREF (2021) as class C0, being the most stable and reliable stations in the EPN. The other four stations are class C1 and are still very stable. The CORS set was completed with the MAS1 and NOT1 stations in the IGS14 Core network (International GNSS Service), resulting in fifteen reference stations (Fig. 1).

Data processing strategy
GPS data are processed using a three-step approach similar to the approach described by Dong et al. (1998) and Herring et al. (2018). In the first step, daily loosely constrained estimates of station coordinates and their covariances are obtained from doubly differenced GPS phase observables (GAMIT) or from PPP using undifferenced observables from just a single GPS receiver (GipsyX). In the second step, time series (glred) in a consistent reference frame (IGb14) are generated by minimizing each day for the selected CORS set (glorg). Outliers and offsets due to earthquakes or equipment changes, as well as periodic signals, are corrected to obtain clean time series. Additionally, in this step, the loosely constrained solutions from GAMIT and GipsyX were combined for each day, generating their time series. In the third step, all data are combined (globk) into a single solution applying generalized constraints (glorg) in the position and velocity to define a uniform reference frame while estimating an overall network translation and rotation (Helmert parameters). Both the second and third steps follow an iterative process. The first time series was stabilized (reference frame definition) using only the well-defined CORS, but the time series that best represents the final velocity solution is the one that is stabilized using all sites in the solution (Herring et al. 2015, 2018) (Figs. S9-S21 in S1_1). Thus, the first velocity solution obtained is used as a priori information for the next iteration, and all sites were used for stabilization. All the time series and velocity solutions are realized in a no-net-rotation (NNR) International Terrestrial Reference Frame (IGb14). Finally, positions, velocities, and uncertainties with respect to the Eurasian-and Nubian-fixed reference frames were estimated by fixing the set of CORS located on the stable part of the plates (glorg). Figure 4 depicts the process flow diagram of the processing strategy described above. Some particular aspects of each strategy are discussed below. Details of the GAMIT and GipsyX processing specifications are presented in Tables S2 and S3 (SI_1).

GAMIT strategy
The Topo-Iberia stations (25) plus the additional CORS (15) are processed simultaneously in only one regional network (40) using GAMIT version 10.71 (Herring et al. 2018). Daily station positions and atmospheric delays are estimated from an ionosphere-free linear combination of GPS phase observables using the double-differencing technique to eliminate phase biases in the satellite and receiver clock oscillators. Some of the main models and processing parameters selected in this strategy are shown in Table 1. The GAMIT solutions are obtained in daily ASCII h-file format that are passed to GLOBK, where the final solution is determined. The integer phase ambiguities are successfully resolved (fixed) with an average of 96.1% for wide-lane (WL) and 92.8% for narrow-lane (NL) and a normalized root-meansquare (NRMS) average of 0.18, which is an indicator of the quality of the solutions. A more detailed explanation of the resolved ambiguities is provided in SI_1 (Fig. S2).

GipsyX strategy
Each of the 40 sites is separately processed with GipsyX version 1.7 (Bertiger et al. 2020) in PPP mode using ionosphere-free carrier phase and pseudorange observables. Some models and processing parameters are common to those used in GAMIT (Table 1). In this case, the integer phase ambiguities are successfully resolved with an average of 99.8% for all RINEX files, with 3832 the average number of the double differences that pass certain established criteria for ambiguity resolution. The GipsyX solutions are obtained in daily gdcov-file format, which is translated to the Solution Independent Exchange (SINEX) format file (IERS 2006) before being passed to GLOBK. A program called gdcov2sinex.sh was specifically created for this conversion. To make these SINEX files loosely constrained, the methodology described by Herring et al. (2015Herring et al. ( , 2016 is followed. Details about the ambiguity resolution and the loosening process are described in SI_1 (Fig. S3). The conversion process and program code are provided in Supplement 2 (SI_2). For ease of use, the executable program file (.sh) is provided in Supplement 3 (zipped folder SI_3.zip).

Combined strategy
As described above, different models and processing parameters are used in GAMIT and GipsyX. Consequently, the noise models from both analyses differ, resulting in different position estimate standard deviations. For the individual final solutions to be comparable, their covariance matrices are rescaled to weight the solutions equally following the combination procedure established by Herring et al. (2016). Two modifications with respect to the original methodology

Considerations for the vertical component
The precision of the GNSS vertical component is generally approximately 2-3 times lower than that of the horizontal component for various reasons, including the geometry and distribution of satellites and phenomena such as ocean loading and atmospheric delays. To improve the accuracy, some changes have been made to the default GAMIT configuration. The default mapping function used for modeling and estimation of the atmospheric delay is the global mapping function (GMF) from fitting numerical weather model (NWM) data over 20 years, sufficient for all but the most accurate studies of heights and for meteorological studies (Herring et al. 2018). However, the VMF1 model is used, providing a more accurate reconstruction of the NWM data. Since having the most accurate value of pressure, from local surface pressure measurements, for the a priori zenith hydrostatic delay (ZHD) was not possible, the next most accurate source was chosen, the VMF1 grid. In addition, zenith delays were estimated every 2 h, although there is little increase in the accuracy of station heights, which greatly increases the processing time. Regarding ocean tide loading, the model with the best possible accuracy at the date of the study was used, the FES2014b model. The default GipsyX configuration was also changed with regard to the VMF1 and FES2014b models.

Time series analysis
The loose solutions previously obtained are considered individually using glred to generate the position time series of each station, with day-to-day repeatability, throughout the span considered in an adequate framework. The raw time series are analyzed to estimate and remove outliers, periodic terms, and discontinuities in an iterative process described in SI_1. The analysis is realized using the tsfit GLOBK program, a command-line tool for fitting time series and generating statistics, and the tsview interactive MATLAB-based program (Herring 2003). Files generated with tsfit/tsview containing outliers, discontinuities, and periodic signals are input to GLOBK to obtain a combined solution and extract more refined positions and velocities. Finally, the noise in the GPS position time series is assumed to be a combination of white noise plus random-walk (RW) noise, which is accounted for using the realistic-sigma algorithm or First-Order Gauss-Markov Extrapolation (FOGMEX) algorithm (Herring et al. 2015;Floyd and Herring 2020). Raw and cleaned time series are shown in Figs. S9-S21 (SI_1).

Analysis of the final results and discussion
Once the GPS data processing described above was completed, estimates of the velocity field based on the final position time series analysis of the data were generated. In the analysis of the results, a distinction was made between the final position time series and the velocity solutions obtained with each strategy, and a comparison between them was made. It is important to note that the FOGMEX algorithm has not yet been applied to these final position time series, which are those stabilized with all the sites in the solution, i.e., the time series that best represent the final velocity solution. The combined solution corresponds to Figs. S9-S21 in S1_1. Finally, the new 3D velocity field of the Topo-Iberia network is compared with previous publications, and the main differences and similarities are highlighted.

Comparison of time series
The statistics (mean and median) of the WRMS position residual scatter in the final fitting of the time series are similar in the three strategies, and the differences between them are on the order of a tenth of a mm. The means of the WRMS scatters of position in the north and east components in the time series for the three strategies are in the range of 0.6-0.7 mm and of 2.0-2.1 mm in the up component, indicating good repeatability. Comparison details are shown in SI_1 ( Fig. S22 and Table S6).

Comparison of velocity solutions
Uncertainties for rates quoted in the text are one-sigma and those shown in the figures are 95% confidence intervals unless stated otherwise. A previous comparison between GAMIT and GipsyX solutions without applying the realistic-sigma algorithm was used to check their degree of matching. Details about this comparison and the different htoglb options to loosen the constraints of the GAMIT and GipsyX solutions are described in SI_1 (Fig. S23). A comparison of the mean and RMS of the differences between the final velocity estimates obtained using the three different strategies is shown in Table 2. The FOG-MEX algorithm is applied to these final velocity solutions to account for the long-period temporally correlated noise and obtain more realistic uncertainties. The agreement between the velocity estimates from each strategy is very good overall, with RMS of differences < 0.2 mm/yr in the horizontal component and < 0.3 mm/yr in the height component. The comparison of velocity estimates from the GAMIT and GipsyX solutions shows that they match at the level of 0.2 mm/yr in the north and east components and 0.3 mm/yr in height. These values are very similar to those resulting from comparing the Combined and GipsyX solutions, but when the Combined and GAMIT solutions are compared, the RMS of the differences decreases to 0.1 mm/yr in the horizontal and up components.
The final horizontal velocity solutions in a Eurasianfixed reference frame are shown in Fig. 5. The Combined solution tends to have an intermediate position between the other velocity solutions, although closer to the GAMIT solution, as shown in the statistics in Table 2. The mean horizontal velocity uncertainty is 0.1 mm/yr for the three strategies. The Combined horizontal velocity solution is presented as the new horizontal CGPS velocity field of the Topo-Iberia network based on 11 years (2008-2019). The velocities obtained for all stations used in this analysis are shown in Tables S7-S9 (SI_1). In addition, to facilitate plotting these velocities with GMT, an ASCII text file (hz_vel.dat) is provided in the zipped folder SI_3.zip. The expected velocities for the reference stations are known from data published by international organizations, such as the IGS and EPN. These velocities and coordinates have been used as a priori values so that the adjustment level indicates the validity of the results. The RMS of the adjustment in the horizontal velocity is < 0.2 mm/yr, denoting an excellent fit (see details in SI_1, Figs. S26-S28).
The final vertical velocities depicted in Fig. 6 confirm that the values from the GAMIT and Combined strategies are quite similar, as indicated by the RMS of their differences (0.1 mm/yr in Table 2). In contrast, in the GipsyX strategy values, there are more differences in magnitude than the other strategies and the opposite results in 25% of cases. The vertical velocities are generally low, with a mean absolute value of 0.5 mm/yr and a mean uncertainty of 0.2 mm/yr for all strategies. The uncertainties are shown in Fig. 7, and a decreasing trend is observed as the processed data span increases in a similar way for the three strategies, coinciding with their linear regression lines. The highest uncertainty in the graph corresponds to the TIOU station, with a data span of 8.5 years. Fifty percent of the Topo-Iberia stations have a data span of over 9.5 years, and according to the graph, this corresponds to an uncertainty of ~ 0.2 mm/yr, which is quite acceptable considering the low vertical velocities estimated.
The Combined vertical velocity solution is shown in Fig. 8 as the first vertical CGPS velocity field of the Topo-Iberia network. The vertical velocity solutions realized in the different reference frames are similar. Generalized subsidence is observed, except in a few cases, with the BENI and LNDA stations showing uplifts of 0.6 ± 0.5 and 0.6 ± 0.4 mm/yr, respectively. Coincidentally, the processed data spans from both stations are the shortest, 4.5 years (Fig. 6), and possibly a longer data span would be needed to obtain lower uncertainty and verify the uplift. Furthermore, the BENI station has been dismantled, resulting in a greater lack of information regarding what happened in its surroundings; hence, although the data recorded at this station were used in this study, the velocities from the BENI station cannot be presented herein with total confidence. Garate et al. (2015) realized the velocity solution in the Eurasia fixed ITRF2008 frame; however, the new velocity solution was realized in IGb14 with respect to Eurasianand Nubian-fixed reference frames. To compare the different fields, the ITRF2008 velocity solution is rotated and translated into the reference frame of the new velocity solution using the velrot v.1.01 program. The resulting average RMS of the fit based on 33 common stations is 0.4 mm/yr, indicating a good adjustment. The new velocity solution is also compared with the velocity field presented by Koulali et al. (2011), although in this case, there are no Topo-Iberia stations in common. In this case, the RMS fit based on 11 common stations is 0.2 mm/yr, also indicating a good adjustment.

Comparison with previous publications
The new horizontal CGPS velocity solution is in general in good agreement with velocities published by Garate et al. (2015). The smallest residual displacements with respect to Eurasia are observed in northern Iberia. Their existing displacement values range from 0.1 mm/yr at ASIN in the northeast to 1.1 mm/yr at EPCU in the central zone (Table S8 of SI_1 and Fig. 9). The stations with the highest displacements are those deployed in Morocco and located on the Nubia plate, showing general northwestward displacements ranging from 3.7 to 4.9 mm/yr. The main differences are visible in the Betic Cordillera, specifically at the ALJI, LIJA, LOJA, and NEVA stations (Fig. 10 top), located in the Betic-Rif deformation zone.
Most of the displacements in the central and eastern Betics are westward, with a variable southwestward component. The AREZ and TGIL stations maintain the westward and southwestward displacements, respectively, presented in previous studies, but their velocities are even reduced by 76% in the case of TGIL, both being below 0.9 mm/yr. The PILA station (1.4 mm/yr) shows northwestward displacement in the eastern Betics, where part of the shortening between the Nubia and Eurasia plates is accommodated (Borque et al. 2019). The NEVA station (1.7 mm/yr) continues showing a displacement roughly parallel to the nearby GRAN station and the farthest CAAL station, as in Galindo-Zaldivar et al. (2015), coinciding exactly with the azimuth of the PALM station (2.8 mm/yr), but the stations as a whole show a higher southwestward component instead of purely westward motion. The LOJA station (2.5 mm/yr), in addition to presenting a reduction in velocity similar to the NEVA and PALM stations, also presents a higher southwestward component integrated into both stations and together with the CABR and MOTR stations of the velocity field of Koulali et al. (2011).
As previously published, displacements in the western Betic Cordillera show westward motion (Garate et al. 2015;Gonzalez-Castillo et al. 2015). At the LIJA station (2.5 mm/ yr), there are no evident changes in motion due to its current uncertainty showing a westward displacement similar to the ROND and LEBR stations. Between the displacements of the Topo-Iberia stations in the western Betic, a maximum is reached at the ALJI station (3 mm/yr), showing present-day west-northwest ward motion rather than westward motion parallel to the LIJA station. The new displacement at the ALJI station, compared with the velocity field of Koulali   2011), is in harmony with the ALGC and UCAD stations, being even more evident with respect to the Nubianfixed reference frame (Fig. 10 bottom). This new direction supports the conclusion of Gonzalez-Castillo et al. (2015) that a clockwise block rotation occurs at present in the northern branch of the tectonic belt Gibraltar Arc. According to the new data presented in this study, this region is currently undergoing a northwest-southeast convergence rate of 3.5 ± 0.1 mm/yr with respect to Nubia. The displacements as a whole are also in agreement with the tectonic model involving rollback and indentation tectonics in the Gibraltar Arc proposed by Galindo-Zaldivar et al. (2022).
It might be believed that the origin of the discrepancies described above could be related to different criteria followed by the authors in the time series analyses, but the case of the NEVA station is remarkable, presenting a new direction of displacement. The NEVA time series have been very stable with the different criteria used, keeping the result constant, without information gaps in the metadata and without unknown strange behaviors during the long activity period (10.5 years). Therefore, the new southwest displacement at the LOJA station, more or less parallel to the NEVA and PALM stations, can be considered coherent. In addition, the data span of the LOJA and ALJI stations analyzed in this work has been extended from 4 to 9.5 years, allowing a more precise determination of their present-day displacements.
With respect to the vertical velocity field, in addition to the higher uncertainties associated with this component, it must be considered that the results presented by different publications usually come from methodologies differing in, among others, processing software, noise Fig. 8 Map of vertical velocities of the Topo-Iberia CGPS network and CORS used to define the framework (names ending with an asterisk). Positive and negative values correspond to uplift and subsidence, respectively. One-sigma uncertainties are represented by vertical error bars model or stations defining the reference system. As a consequence, the vertical velocity estimates presented tend to have quite a few differences. All of this makes a detailed comparison of the velocity at each station less conclusive than whether or not there is an agreement between regions of subsidence or uplift. On the other hand, the density of the Topo-Iberia network and CORS used in this study, although covering the study area, does not allow for establishing rates of subsidence or uplift with confidence for the entire Iberian Peninsula and Morocco. It should also be noted that the vertical velocities determined are so small that they may not imply subsidence or uplift. However, as in previous publications, general subsidence is observed in the Iberian Peninsula (mean value of 0.5 ± 0.1 mm/yr) except at some isolated stations (Fig. 11). A mean subsidence velocity of 0.6 ± 0.1 mm/yr is located in the southern Iberian Peninsula and approaches zero close to the southern Pyrenees, in agreement with Serpelloni et al. (2013), although they establish subsidence rates ~ 2 mm/yr higher than those in this work. They did not use Topo-Iberia stations, and although they present a dense spatial distribution of stations, they only covered the central and southeastern regions. The velocity rates close to zero in the Pyrenees also coincide with those presented by Nguyen et al. (2016), where there are four stations in common and the estimated velocity coincides at two of them (LNDA and FUEN), which could confirm the uplift at the LNDA station. Fig. 9 Comparison between the new horizontal CGPS velocity field (95% confidence interval) and previous publications with respect to a Eurasian-fixed reference frame. Note that a factor of 0.5 has been applied to the velocity field uncertainties of Koulali et al. (2011) for figure clarity. The dashed black box represents the region enlarged in Fig. 10. Similar comparison with respect to Nubia is shown in Fig. S25 (SI_1) 1 3 155 Page 12 of 15

Summary and conclusions
Topo-Iberia network data available up to 2019, with 50% of stations having over 9.5 years of data, have been processed, presenting a new 3D velocity field based on 11 years.
A quality check was carried out on the entire network, from which the optimal data spans to be processed were determined. In addition, a global vision of the data quality from the Topo-Iberia stations is presented for the first time. Fig. 10 Comparison between the new horizontal CGPS velocity field (95% confidence interval) and previous publications with respect to Eurasian-(top) and Nubian-fixed (bottom) reference frames. Note that a factor of 0.5 has been applied to the velocity field uncertainties of Koulali et al. (2011) Page 13 of 15 155 GAMIT/GLOBK and GipsyX were used with the intention of having different analyses to compare and avoid possible errors introduced during processing while maintaining the premise of using different GNSS processing techniques. To make the individual solutions comparable, the methodology described by Herring et al. (2016) was followed, but in this case, only the GipsyX solution, and not the GAMIT solution, was rescaled by a factor obtained from the relation between the sigma values in the time series of both solutions. The final velocity solutions were not combined or averaged, but the two independent loosely constrained solutions from GAMIT and GipsyX were merged using a kalman filter (GLOBK). This required converting the GipsyX solutions to SINEX files, and a new program specifically created for it called gdcov2sinex.sh is provided.
The means of the WRMS scatters of position in the north and east components in the time series for the three strategies are in the range of 0.6-0.7 mm and of 2.0-2.1 mm in the up component, indicating good repeatability.
A new horizontal CGPS velocity field of the Topo-Iberia network based on 11 years (2008-2019) has been presented. The mean horizontal velocity uncertainty (1σ) is low, 0.1 mm/yr. The new velocity field is in good agreement with the results of previous publications, except for some discrepancies in the Betic Cordillera, where the new horizontal velocities are lower. The smallest and highest residual displacements with respect to Eurasia occur in northern Iberia and Morocco, respectively. The stations in the Betic central sector present a higher southwestward component instead of purely westward motion. The main discrepancy in the western Betic Cordillera is observed at the ALJI station, where  Table S2, and Nguyen et al. 2016 Table 3). A black triangle symbol in the eastern Iberian Peninsula exceeds the color scale bounds and corresponds to a velocity value of −9.3 mm/year a maximum is reached (3 ± 0.1 mm/yr), showing presentday west-northwestward rather than westward motion. The results support the conclusion of Gonzalez-Castillo et al. (2015) that a clockwise block rotation occurs at present in the northern branch of the Gibraltar Arc. The velocities indicate a northwest-southeast convergence rate of 3.5 ± 0.1 mm/ yr with respect to Nubia. The displacements as a whole are also in agreement with the tectonic model involving rollback and indentation proposed by Galindo-Zaldivar et al. (2022).
A vertical CGPS velocity field of the Topo-Iberia network has been presented for the first time. The vertical velocities are low for the three strategies (a mean absolute value of 0.5 mm/yr and a mean uncertainty of 0.2 mm/yr). Considering the limitations of the vertical component estimation and that the vertical velocities determined are so small that they may not imply subsidence or uplift, general low subsidence is observed in the Iberian Peninsula, as in previous publications, except at some isolated stations, such as the LNDA station, which shows an uplift of 0.6 ± 0.4 mm/yr. A mean subsidence velocity of 0.6 ± 0.1 mm/yr is located in the southern Iberian Peninsula and approaches zero close to the southern Pyrenees. The results offered additional information on the vertical velocity in the Iberian Peninsula and Morocco, complementing the few studies that have estimated it.