*Article* **Quantitative Evaluation of Environmental Loading Products and Thermal Expansion Effect for Correcting GNSS Vertical Coordinate Time Series in Taiwan**

**Bin Liu 1,2, Xiaojun Ma 1,2, Xuemin Xing 1,2,\*, Jianbo Tan 1,2, Wei Peng 1,2 and Liqun Zhang 1,2**


**Abstract:** We explored the driving factors of nonlinear signals in vertical coordinate sequences of stations in a Taiwan global navigation satellite system (GNSS) network, including atmospheric loading (ATML), hydrological loading (HYDL), and non-tidal ocean loading (NTOL) effects. At the same time, we used the finite element analysis software MIDAS to quantify the vertical displacements of different types of monuments due to the thermal expansion effect, including deep drilled braced (DDB) and short drilled braced (SDB). By quantitatively comparing the correction results of GNSS time series with different single mass loading models, we found that there was little difference in the correction of different environmental loading products. We compared different combinations of each loading product to correct the GNSS time series, and finally selected the best combination suitable for Taiwan GNSS network, that is, ATML (GFZ\_ECMWF IB) + HYDL (IMLS\_MERRA2) + NTOL (IMLS\_MPIOM06). We found that the spatial and temporal models of ATML and NTOL are very similar, with non-tidal atmospheric loading and non-tidal ocean loading working together, a pattern that may be related to tropical cyclones. Both models also showed good correction effect on GNSS stations in the western plain of Taiwan, but with limited correction effect in the eastern part of Taiwan. This may be due to the influence of the subtropical monsoon climate in Taiwan and the barrier of the central mountain range, resulting in obvious differences between eastern and western Taiwan. The hydrological loading was found to act in the opposite way to the thermal expansion effect in the temporal domain, indicating that some displacements in hydrological loading may cancel out displacements caused by the thermal expansion effect. This aspect of displacement is not included in the hydrological loading model but should be considered when accurately estimating the temporal and spatial variation of water storage capacity in Taiwan using GNSS observed displacements.

**Keywords:** GNSS; environmental loading; thermal expansion effect; monument types

#### **1. Introduction**

Moving masses of surface fluid layers interact with the solid Earth, resulting in periodic weak deformation and surface displacements [1]. Permanent reference stations on the ground can capture these non-tectonic deformation signals by long-term, continuous, and high-precision monitoring. Previous studies have identified three main classes of nonlinear changes causing the nonlinear variation in Global Navigation Satellite System (GNSS) derived coordinate time series. The first is error caused by imperfect GNSS data processing models or other systematic errors, such as high-order ionospheric delay errors, multipath effects [2–5], etc. The second is due to real nonlinear movement of the reference stations, including periodic displacements of the reference stations caused by geophysical effects, such as environmental loading (atmospheric loading, hydrological loading, or nontidal ocean loading) and thermal expansion effects of bedrock and the monument [6–11],

**Citation:** Liu, B.; Ma, X.; Xing, X.; Tan, J.; Peng, W.; Zhang, L. Quantitative Evaluation of Environmental Loading Products and Thermal Expansion Effect for Correcting GNSS Vertical Coordinate Time Series in Taiwan. *Remote Sens.* **2022**, *14*, 4480. https://doi.org/ 10.3390/rs14184480

Academic Editors: Alex Hay-Man Ng, Linlin Ge, Hsing-Chung Chang and Zheyuan Du

Received: 8 August 2022 Accepted: 5 September 2022 Published: 8 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

exhibited as seasonal periodic oscillation in the GNSS coordinate time series. The third is the periodic changes caused by observation noise. By studying the periodic deformation of environmental loads, the annual periodic amplitude of non-tectonic deformation of GNSS can be effectively reduced for improving the velocity estimation accuracy.

There has been increased attention paid to quantitative analysis of the contributions of different factors to GNSS coordinate time series, especially for the environmental loading effect. However, the correction effect of environmental loading deformations for GNSS coordinate time series varies greatly in different regions [11–14]. Moreover, for the quantitative calculation of loading deformation and displacements, the use of different Earth models, reference frames, mass loading data models, and coastline data from different institutions also give different calculation results. Li et al. [13] studied the individual contributions of five different continental water storage loading models, six different non-tidal atmospheric loading models, and five different non-tidal ocean loading models to the 220 GNSS coordinate time series in the Crustal Movement Observation Network of China (CMONOC), considering the joint contributions of different permutations and combinations of load products to GNSS time series. The results have a RMS reduction in the vertical component of 20%. Wu et al. [14] compared the contributions of environmental loading sequences from the German Research Center for Geosciences (GFZ) and the School and Observatory of Earth Sciences (EOST) to GNSS vertical time series. Compared with the RMS reduction value for individual stations, the combined model using sequences from both GFZ and EOST was better than use of GFZ alone to correct coordinate time series. Results can also be influenced by the independent climatic and topographic characteristics of different regions. Jiang et al. [11] considered the differences of hydrological data, grid interpolation, and terrain correction in the optimum model data (OMD) model for reduced scattering at 74% of global stations. One suggestion from the above results is that the loading sequences from different organizations should be adapted to local conditions, and the correction of GNSS stations coordinate time series should be carefully selected for different regions.

In addition to environmental loading, the thermoelastic deformation of the antenna monument of GNSS reference station and its bedrock driven by periodic temperature changes also contribute to nonlinear signals in time series data [7,15]. On the seasonal scale, the annual amplitude caused by the thermal expansion effect of the monument can exceed 5 mm for GNSS reference stations where the monuments are located at high latitudes [10]. Dong et al. [7] found that the annual vertical deformation of the surface caused by a change of surface temperature is less than 0.56 mm based on the semi-infinite space model. Yan et al. [16] calculated the influence of temperature change on the underground bedrock and cement pier of GNSS station in China, and showed that the maximum annual amplitude of vertical displacements caused by thermal expansion effect can reach 2.8 mm. Based on the three-dimensional full space thermoelastic deformation model, Tan et al. [17] calculated the influence of temperature change on the three-dimensional annual amplitude for China. The results showed that temperature change was an important factor causing continental surface deformation in addition to mass loading. Jia et al. [18] quantitatively analyzed the difference of vertical deformation between GNSS and Gravity Recovery and Climate Experiment (GRACE), and analyzed the influence of thermal expansion effect on vertical displacements of GNSS and estimation of vertical loading deformation of surface by GRACE. The results showed more consistent vertical deformation of GNSS and GRACE for the Chinese mainland after correcting for the thermal expansion effect in GNSS time series data. Boccara et al. [19] compared the differences between numerical models of land surface temperature from different sources at the global scale, and found that the annual amplitude of land surface temperature data provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) was closest to the measured temperature data. Vertical displacements due to thermal expansion effect include the monument (TEM) on the ground and the bedrock (TEB) underground, but not any displacements caused by different specifications and materials of the monument. The GNSS network in Taiwan Province, China is operated by multiple organizations, so Rao et al. [20] classified more

than 90% of the base types of stations and studied the performance of GNSS coordinate time series for the different bases. They found that the coordinate time series of metal rod monument has significant annual periodic variation, and the months of annual amplitude uplift and subsidence are similar to temperature variation.

Taiwan Province is located in an active tectonic area. With continuous collision between Luzon arc and Chinese mainland margin, there are many thrust faults and folds. Since 1990, there have been many destructive earthquakes with magnitude greater than 6. To obtain a more reliable velocity field, there is a significant need for more accurate estimation of nonlinear signals such as environmental loads and geophysical factors in GNSS coordinate sequence data. Kumar et al. [21] studied the sources of GNSS commonmode errors (CME) in Taiwan Province, and showed that 90% of the CME variations in Taiwan are related to non-seasonal ATML displacements. Ma et al. [22] used independent component analysis (ICA) to separate the signals related to atmospheric and hydrological loading from GNSS residual coordinate time series in Taiwan Province. Lai et al. [23] used data from GNSS stations in Taiwan to study inversion terrestrial water storage (TWS). The results showed that GNSS shows large annual vertical displacements in the plain area of southwest Taiwan and the Huadong Valley in the east of Taiwan. Using GNSS, GRACE, and hydrological (precipitation, GLDAS and LSDM assimilation models, and in situ groundwater level) datasets, Hsu et al. [24] systematically studied the spatial and temporal changes of water storage capacity in Taiwan. The results showed basically consistent spatial patterns of annual water storage capacity estimated by GNSS, GLDAS, and precipitation data, indicating that there is significant seasonal hydrological loading fluctuation in Taiwan.

Based on the findings of these previous studies, we analyzed the influence of atmospheric loading (ATML), hydrological loading (HYDL), and non-tidal ocean loading (NTOL) on Taiwan GNSS time series data. We considered the influence of temperature change on different bases of the stations, and accurately calculated the change in vertical displacements due to temperature changes of GNSS stations. We studied 84 GNSS sites in Taiwan Province from 1 January 2008 to 31 December 2012. These data were corrected using six ATML, five HYDL, and three NTOL products provided by International Mass Loading Service (IMLS), GFZ, and EOST. Displacements caused by thermal expansion effect (TEE) were characterized. Finally, the RMS reduction value was used to evaluate the influence of non-tectonic signals on GNSS sequences. By comparing the mass loading models provided by different institutions, we determined a group of combined mass loading models suitable for correcting GNSS time series in Taiwan Province. At the same time, we quantified the vertical displacements of DDB and SDB monument stations affected by temperature changes. These would be helpful to study the driving factors in nonlinear signals in GNSS vertical coordinate time series in Taiwan Province.

#### **2. Data and Methodology**

#### *2.1. GNSS Data*

The Taiwan continuous GNSS Network was established by the Institute of Earth Sciences (IESAS) in 1990. Currently, it consists of 539 stations that belong to different institutions and are operated by the Central Weather Bureau (CWB), IESAS, the Central Geological Survey (CGS), and the Ministry of the Interior (MOI). We utilized the classification of stations base forms devised by Rao et al. [20] and classified the stations into six types (Figure 1): deep drilled braced (DDB) monument (38), short drilled braced (SDB) monument (217), concrete pillar (30), roof (188), metal rod (14), and other forms (52).

After the above classification of station bases, we selected 84 DDB and SDB stations with high data quality for follow-up analysis. DDB and SDB stations were selected because these were constructed according to fixed specifications of the Central Weather Bureau and the Southern California Integrated GPS Network (SCIGN) (Figure 2a,b). DDB and SDB stations differ in that the DDB base contains five 12 m 316 stainless steel pipe supports that are buried 10 m deep underground. SDB contains only four stainless steel pipe supports, buried 1 to 2 m underground. Due to limited funds, regional restrictions, and construction difficulties in high mountain areas, GNSS stations in Taiwan have mostly adopted the SDB pedestal. The advantage of this pedestal is that these stations are mainly distributed in suburbs and mountainous areas far away from cities, and this kind of pedestal is relatively stable and is less noisy. Although roof pedestals account for 35% of GNSS observation network in Taiwan, these stations are mostly located in cities and areas with many buildings. To calculate the thermal expansion effect, many factors must be considered, so roof stations were not considered in this study.

**Figure 1.** Distribution map of GNSS stations monument types in Taiwan Province. The circle represents the stations that are still in operation, the pentagram represents the stations that have been deactivated, and the triangle represents the stations used in this paper. Different colors represent different monument types: red are the deep drilled braced (DDB) monument stations, blue are the short drilled braced (SDB) monument stations, cyan are the concrete pillar monument stations, magenta are the roof stations, and purple are the metal rod type stations.

The GNSS data from 2008 to 2012 is processed using the GAMIT/GLOBK software packages [25], and a composite daily solution in the ITRF 2008 reference framework [26] is generated for subsequent time series analysis. IERS 2003 model and FES 2004 models were used for data processing to correct for polar tide and sea-ocean tidal loading [27,28]. Separate documents for each month provided by different institutions. After cubic spline interpolation, we obtained the predicted time series of loading products at each GNSS station. Our goal was to study the influence of environmental loading and thermal expansion on GNSS time series. After eliminating outliers, we used the following standard procedure to fit the time series data [29–31], and remove the trend term and offset terms of the coordinate time series.

$$y(t\_i) = a + bt\_i + \varepsilon \sin(2\pi t\_i) + d\cos(2\pi t\_i) + \varepsilon \sin(4\pi t\_i) + f\cos(4\pi t\_i) + \sum\_{j=1}^{n\_\mathcal{S}} g\_i H(t\_i - T\_{\bar{\mathcal{g}}i}) + \varepsilon\_i \tag{1}$$

where *ti*(*i* = 1, 2, ··· , *N*) is the daily solution epoch units of years, *a* and *b* are the intercept and linear rate, *c* and *d* are annual term coefficients, *e* and *f* are semi-annual coefficients, *H ti* <sup>−</sup> *Tgi* is a Heaviside step function, *ε<sup>i</sup>* is the model residual.

**Figure 2.** Design drawing of instrument monument. (**a**) It is a deep drilled braced monument station. The DDB monument is composed of five pillars, each of which extends 35 feet underground and 62 inches on the ground. (**b**) It is a short drilled braced monument station. The SDB monument is composed of four pillars, each of which extends 60–72 inches underground and 40–55 inches on the ground. Refer to UNAVCO website for more details (http://www.unavco.org/ (accessed on 7 August 2022)).

#### *2.2. Surface Mass Loading Effects*

Independently derived global models of major seasonal mass loading processes were used to assess the specific contribution of these load sources to the vertical data of Taiwan Province's GNSS observation network from 2008 to 2012. In this study, we focused on the seasonal movements in the vertical direction. We used the reanalysis numerical weather model MERRA2 and the operational semi-frozen numerical weather model GEOSF-PIT provided by the International Mass Loading Service (IMLS, http://massloading.net/ #Download (accessed on 7 August 2022)) [32–35] to calculate the atmospheric mass loading displacements and the hydrological loading displacements, respectively. Non-tidal ocean loading displacements were calculated using the MPIOM06 model provided by the Max Planck Institute for Meteorology and GFZ. The load products provided by IMLS are built on a regular 2 × 2 global grid, and pressure loading data corresponding to any point on Earth are available online on demand. The School and Observatory of Earth Sciences (EOST, http://loading.u-strasbg.fr/ (accessed on 7 August 2022)) [36] provides global displacement maps with a resolution of 0.5 degrees, with predicted displacements provided in ascii and netcdf data formats in the center of figure (CF) and center of mass (CM) reference frames. The ATML were loaded by modified inverted barometer model (ATMIB), induced oceanic loading (ATMMO), and hydrological loading including GLDAS/Noah models, with ERA interim 6-hour sampling by adjusting monthly model and deducting air tides (ERAin). Non-tidal ocean loading includes estimates of the circulation and climate of the ocean (ECCO1) and ECCO2 (follow-on ECCO, phase II). The German Research

Center for Geosciences (GFZ, http://rz-vm115.gfz-potsdam.de:8080/repository (accessed on 7 August 2022)) [37] provides non-tidal atmospheric loading products for atmospheric surface pressure enforced by ECMWF, hydrological loading products for calculations using the Land Surface Discharge Model (LSDM), and non-tidal ocean loading products for the EMPIOM ocean circulation model for atmospheric data operated by ECMWF. GFZ can also provide the calculation of loading products at GNSS sites. The above environmental loading time series were derived in the CF reference frame, as detailed in Table 1.


**Table 1.** Environmental loading product provided by EOST, GFZ, and IMLS.

#### *2.3. Temperature and Precipitation Data*

The temperature data we used were from the reanalyzed datasets provided by the European Center for Medium-Range Weather Forecasts (ECMWF). ECMWF data provide a land surface temperature grid with the highest resolution of 0.125◦ × 0.125◦ , and time resolution of four points per day, measuring 0, 600, 1200, and 1800 h (UTC) (https://apps. ecmwf.int/datasets/interim-full-daily (accessed on 7 August 2022)). In this experiment, we selected 0.5◦ × 0.5◦ resolution temperature grid data at 2 m on the ground at 600 (UTC), sampled every day from 1 January 2008 to 31 December 2012.

The spatial resolution was 0.1◦ × 0.1◦, and the rainfall data sampled month by month were provided by the global precision measurement mission (GPM, https://disc.gsfc.nasa. gov/datasets?keywords=GPM&page=1 (accessed on 7 August 2022)).

#### *2.4. Methodology*

2.4.1. Calculation of Thermal Expansion Displacement

The effect of thermal expansion caused by temperature change on the vertical displacements of GNSS stations was divided into aboveground and underground contributions [10]. Generally, for the length change caused by expansion with heat and contraction with cold of the monument with GNSS antenna above the ground surface, the calculation equation is as follows [10,38]:

$$
\Delta L = \alpha L \Delta T \tag{2}
$$

where Δ*L* is the length change of the monument, the linear thermal expansion coefficient *<sup>α</sup>* = 1.2 × <sup>10</sup>−5/◦C, <sup>Δ</sup>*<sup>T</sup>* is the temperature change, as the difference between the average value in the relative observation period, and *L* is the height of the monument.

DDB and SDB are composed of five and four stainless steel pipes, respectively. Compared with a single cement column, the connection between the bottom of SDB and DDB base rod and the ground is simulated by fixed support, with the top of the rod is rigidly connected. The whole structure is a high-order statically indeterminate structure, and the

additional internal force generated by temperature change can lead to different displacements of the top of the rod from that of a straight rod. Therefore, it is not appropriate to consider the length of cement pier as a separate parameter in Equation (2), but the displacement change caused by expansion with heat and contraction with cold of stainless steel pipe remains proportional to the temperature difference. We used the finite element analysis software MIDAS to model two different types of pedestal models and added a temperature field. We calculated a total vertical displacement of the two models for a temperature difference of +10 ◦C, and calculated the daily vertical displacements of DDB and SDB pedestals according to Equation (3).

$$
\Delta \mathbf{l} = \frac{\Delta T}{10 \, ^\circ \text{C}} \times d\_{SDB/DDB} \,\tag{3}
$$

where *dSDB*/*DDB* is the vertical displacement calculated for DDB or SDB models in MI-DAS software.

The influence of temperature change on bedrock through heat conduction must also be considered, as well as the resulting change in vertical displacement of GNSS station as follows [7]:

$$
\Delta h = \frac{1+\delta}{1-\delta} \alpha A\_T \sqrt{\frac{\kappa}{\omega}} \cos \left(\omega t - \varphi\_T - \frac{\pi}{4}\right) \tag{4}
$$

where *<sup>δ</sup>* ≈ 0.266, the thermal diffusion coefficients *<sup>κ</sup>* = <sup>1</sup> mm2/s, *AT*, *<sup>ω</sup>*, and *<sup>ϕ</sup><sup>T</sup>* represent the amplitude, angular frequency, and initial phase of annual term of temperature change, respectively. The total displacements caused by thermal expansion effect can be obtained by adding Equations (3) and (4).

#### 2.4.2. Evaluation Metrics

To make a quantitative comparison between the GNSS time series and a mass loading displacement, we defined a misfit function for the two datasets [39]:

$$D = \frac{1}{n} \sqrt{\sum\_{i=1}^{n} \frac{\left(c\_{i2}^{-2} + d\_{i2}^{-2}\right) \left[\left(c\_{i1} - c\_{i2}\right)^{2} + \left(d\_{i1} - d\_{i2}\right)^{2}\right]}{\max\left(c\_{i2}^{-2} + d\_{i2}^{-2}\right)}} \tag{5}$$

where *ci*1, *ci*2, *di*1, and *di*<sup>2</sup> are the amplitudes of the harmonic components of the two compared datasets defined in Equation (1), and *n* is equal to 84. A lower value of *D* indicates a better consistency of the periodic components between the two datasets.

We used RMS reduction value to evaluate the correction effect of products with different environmental loading. The RMS reduction value can be expressed as follows [40]:

$$RMS\_{reduction} = \frac{RMS\_{GNSS} - RMS\_{\text{(GNSS-loading)}}}{RMS\_{GNSS}} \tag{6}$$

where *RMSGNSS* and *RMS*(*GNSS*−*loading*) are the RMS values of the vertical GNSS time series before and after ATML, HYDL, NTOL, and TEE correction. A larger *RMSreduction* value represents a better correction effect.

#### **3. Results**

#### *3.1. Quantitative Assessment of Environmental Loading Effects*

#### 3.1.1. ATML Effects

As shown in Figure 3a, taking the ATML provided by IMLS\_GEOSFPIT as an example, the annual vertical movements caused by atmospheric loading are in the range of 1–1.6 mm, with the peak uplift time around July. The amplitudes and phases distribution of GNSS time series are quite messy, and, except, for a few GNSS stations whose annual vertical movements reach 11 mm, most stations are in the range of 3–6 mm. The phases are bounded by a central mountain range, with a peak uplift time in western Taiwan from May to June

and that in eastern Taiwan from February to April. Thus, atmospheric loading has a positive effect on the vertical time series of GNSS stations in western Taiwan, but a negative effect on the stations in eastern Taiwan.

**Figure 3.** Environmental loading models and GNSS time series amplitudes and phases distribution. The arrow in the figure points to the month with the largest surface displacements uplift. The month rotates clockwise, and the length of the arrow represents the amplitude. According to Equation (1), the amplitude *<sup>A</sup>* <sup>=</sup> <sup>√</sup>*c*<sup>2</sup> <sup>+</sup> *<sup>d</sup>*2, annual phase *<sup>ϕ</sup>* <sup>=</sup> tan−1(*d*/*c*). (**a**) ATML (IMLS\_GEOSFPIT). (**b**) HYDL (IMLS\_MERRA2). (**c**) NTOL (GFZ\_EMPIOM).

We used Equation (5) to quantify the temporal and spatial correlation between six ATML products and vertical GNSS time series. The results are shown in Table 2.


**Table 2.** Consistency between CWS products and GNSS vertical time series.

The results of the RMS value reduction of the vertical GNSS stations coordinate time series are shown in Figure 4 and Table 2, where upward red arrows indicate positive correction and the downward blue arrows indicate negative correction. For visualization purposes, we used the adjustable tension continuous curvature splines [41] to interpolate the RMS reductions for each station displayed as a color profile. The RMS values of the vertical GNSS time series corrected by EOST\_ECMWF IB, EOST\_ECMWF, EOST\_ERA interim, GFZ\_ECMWF IB, IMLS\_GEOSFPIT, and IMLS\_MERRA2 were reduced by −4.77–6.29%, −5.95–8.16%, −4.74–6.25%, −3.54–5.94%, −4.29–5.84%, and −4.25–6.16%, respectively, and the mean RMS reduction values of forward correction results were 2.47%, 3.03%, 2.47%, 2.39%, 2.23%, and 2.45%, respectively. According to the correlation between the six ATML products and GNSS data and the reduction of RMS value after correction, ATML (EOST\_ECMWF) products have the best positive effect for the correction of seasonal changes.

**Figure 4.** The RMS reduction value after GNSS time series correction of six ATML products. The color contours show the interpolated of RMS reduction, and the red (up) arrow indicates the positive response of ATML to GNSS, and the blue (down) arrow indicates the negative response. (**a**–**f**) respectively represent EOST\_ECMWF IB, EOST\_ECMWF, EOST\_ERA interim, GFZ\_ECMWF IB, IMLS\_GEOSFPIT, and IMLS\_MERRA2.

#### 3.1.2. HYDL Effects

As shown in Figure 3b, for HYDL provided by IMLS\_MERRA2 as an example, we can see annual vertical movements caused by hydrological loading of 1.6–2.5 mm for the whole station area, with phases around February. Compared with the western part of Taiwan, the eastern part of Taiwan is closer to the hydrological loading phases. Hydrological loading plays a positive role for most stations in Taiwan Province, so hydrological loading has a significant impact on the change of GNSS time series. We used Equation (5) to quantify the temporal and spatial correlation between five HYDL products and vertical GNSS time series, and the results are shown in Table 3.

**Table 3.** Consistency between HYDL products and GNSS vertical time series.


The RMS reduction results are shown in Figure 5 and Table 3. The RMS reduction ranges of vertical GNSS time series corrected by EOST\_ERA interim, EOST\_GLDAS, GFZ\_LSDM, IMLS\_GEOSFPIT and IMLS\_MERRA2 were calculated as −3.84–6.04%, −2.90–4.56%, −5.59–5.43%, −3.60–7.10%, and −4.17–7.29% , respectively, with mean RMS reduction values of forward correction results of 2.40%, 1.97%, 2.25%, 2.72%, and 2.62%, respectively. The GLDAS/Noah model is mainly composed of soil moisture and excludes the influence of surface water such as lakes and rivers. According to the correlation between five HYDL products and GNSS data and the reduction of RMS value after correction, HYDL (IMLS\_GEOSFPIT) products have the best positive effect to correct seasonal changes.

**Figure 5.** The RMS reduction value after GNSS time series correction of five HYDL products. The red (up) arrow indicates the positive response of ATML to GNSS, and the blue (down) arrow indicates the negative response. (**a**–**e**) respectively represent EOST\_ERA interim, EOST\_GLDAS, GFZ\_LSDM, IMLS\_GEOSFPIT, and IMLS\_MERRA2.

#### 3.1.3. NTOL Effects

As shown in Figure 3c, using NTOL provided by GFZ\_EMPIOM as an example, the annual vertical movements caused by non-tidal ocean loading were calculated as 0.5–1.2 mm, with phases (peak uplift) occurring from June to July. Similar to the atmospheric loading model, according to the phase differences between non-tidal ocean loading and GNSS time vertical series. The correction effect is better for western Taiwan than for eastern Taiwan.

We used Equation (5) to quantify the spatiotemporal correlation between three NTOL products and vertical GNSS time series. The results are shown in Table 4. The RMS reductions of the vertical GNSS time series corrected by EOST\_ECCO2, GFZ\_EMPIOM, and IMLS\_MPIOM06 are shown in Figure 6 and Table 4. The RMS values were decreased by −2.59–2.96%, −2.51–6.41%, and − 2.24–5.14%, respectively, with mean RMS reduction values of forward correction results of 0.95%, 2.00%, and 1.87%, respectively. There was obvious difference between the non-tidal ocean loading provided by GFZ\_EMPIOM and EOST\_ECCO2 to correct the vertical time series of GNSS. NTOL (GFZ\_EMPIOM) products exhibited the best positive effect to correct seasonal changes.


**Table 4.** Consistency between NTOL products and GNSS vertical time series.

**Figure 6.** The RMS reduction value after GNSS time series correction of three NTOL products. (**a**–**c**) respectively represent EOST\_ECCO2, GFZ\_EMPIOM and IMLS\_MPIOM06.

#### *3.2. Quantitative Assessment of Thermal Expansion Effects*

We used finite element analysis software MIDAS to model the pedestal above the ground surface for DDB and SDB GNSS stations, and then added a temperature field. For a temperature difference of +10 ◦C, total displacement in the vertical direction was calculated for the two pedestals as 0.252 mm (DDB, Figure 7a) and 0.219 mm (SDB, Figure 7b). Combined with Equation (3), the daily displacements change of each GNSS station was recorded as TEM. According to Equation (4), the vertical displacement change caused by thermal expansion of bedrock was calculated and is denoted as TEB. The phases of TEB lag 45 d behind those of TEM from Equation (4).

As shown in Figure 7c,d, the annual vertical total displacements caused by TEE were calculated as about 0.3–0.6 mm, with peak uplift occurring around August. At this time, the summer temperature is the highest in the northern hemisphere, but Taiwan Province is located at middle and low latitude, so the annual average temperature difference exhibits little change. In addition, because the phases difference between TEE and GNSS coordinate time series is more than 90◦, the RMS values are negative after deducting the TEE. This is similar to the analysis of Yunnan Province of the Chinese mainland, similarly located at middle and low latitude.

**Figure 7.** (**a**,**b**) MADIS software models DDB and SDB monuments and calculates the surface vertical displacements of the two monument bases under the condition of temperature difference of +10 ◦C. (**c**,**d**) The annual average amplitude and annual phase of the total vertical displacements of DDB and SDB monuments caused by thermal expansion effect.

#### **4. Discussions**

#### *4.1. Optimal Combination of Environmental Loading Models*

We quantitatively evaluated the correction of GNSS coordinate time series by three individual mass-induced displacements. Atmospheric loading and non-tidal ocean loading show similar spatial distribution patterns, with peak uplift phases from June to July and slightly smaller annual amplitudes of non-tidal ocean loading compared to atmospheric loading. This is because Taiwan Province is in a subtropical monsoon climate zone, with many tropical cyclones from July to September. Tropical cyclones cause atmospheric pressure changes and also lead to abnormal migration of seawater quality under the action of wind and atmospheric pressure changes. Non-tidal loading caused by tropical cyclones is the result of the combined action of non-tidal atmospheric loading and non-tidal ocean loading [42]. For the time series correction of GNSS stations, non-tidal atmospheric loading and non-tidal ocean loading are obviously corrected at stations in western Taiwan. This is because this region is mostly plains and GNSS stations are densely distributed, but there are many mountainous areas in the eastern part of Taiwan that block the central mountain range in the middle. Due to the influence of topographic factors, the improvement effect of atmospheric loading and non-tidal ocean loading on eastern stations is not good. Hydrological loading is corrected well at most stations in Taiwan. It is worth noting that, compared

with the western Taiwan where the stations are densely distributed, the maximum subsidence and uplift of GNSS stations along the Huadong Valley in eastern Taiwan occurred in August and February, respectively, which is consistent with the precipitation cycle. After correction for hydrological loading, the reductions of RMS values are equivalent to or even greater than for western stations because the Huadong Rift Valley is formed by the deposition of river channels and alluvial fans and the groundwater layer is thicker.

As we can see, there is no obvious overall difference of different models in the three mass loadings for the correlation between amplitudes and phases and GNSS time series or the decrease of RMS values after correction. This is because the GNSS network in Taiwan Province differs from the global or CMONOC networks and belongs to a small-scale area. In addition to the comparison between individual loaded products, we created all possible combinations of different organization loaded products. We analyzed 90 different combinations in total and they are sorted in descending order according to its mean RMS reduction. Table 5 lists the mean RMS reduction for the first five combinations. Finally, we chose ATML (GFZ\_ECMWF IB) + HYDL (IMLS\_MERRA2) + NTOL (IMLS\_MPIOM06) combination as the best combination mass loading model (CML) suitable for GNSS network in Taiwan. We then calculated the amplitudes and phases of CML according to Equation (1). As shown in Figure 8a, the annual vertical movements caused by CML were calculated as 2–3 mm, and uplift phases were concentrated in mid-May. CML modification of the vertical coordinate time series of GNSS stations is shown in Figure 8b, with the maximum forward RMS reduction of 13.87% (SHLU), and the mean RMS reduction of 5.28%.

**Table 5.** Vertical RMS reduction statistics for selected combinations.


**Figure 8.** (**a**) Average annual amplitude and annual phase distribution of combined mass loading model. (**b**) The positive RMS reduction value of GNSS time series corrected by combined mass loading model.

#### *4.2. The Effect of the Thermal Expansion*

In Section 3.2, we calculated the vertical displacements of SDB and DDB pedestals caused by thermal expansion. The annual average amplitude of vertical displacement is 0.5 mm, accounting for 1.67% of the annual average amplitude of GNSS vertical displacement and about 10% of the annual average amplitude of hydrological loading displacement. Because of the phase difference, the GNSS time series increases after deducting the thermal expansion effect. This may be because the displacements are out of phase, which could be caused by the migration of some surface mass and the thermal expansion effect (Figure 9a), leading to a lack of loading displacements. Lai et al. [23] used GNSS to infer that the terrestrial water storage (TWS) reached the annual maximum and minimum values around August and February, respectively. This is the same timing as the phase mode of vertical displacements caused by thermal expansion effect, but out of phase with the hydrological model. Hsu et al. [24] used geodetic (GNSS and GRACE) and hydrological (precipitation, GLDAS and LSDM assimilation models and groundwater level) datasets to systematically study temporal and spatial water storage changes in Taiwan, and found significant seasonal water load fluctuation in Taiwan, but the total water storage in hydrological assimilation model was underestimated compared with that obtained by GNSS estimation. The difference comes from the complexity of transient water storage due to changes in rainfall patterns, infiltration rates, soil saturation, and runoff.

**Figure 9.** (**a**) Mean hydrological loading displacements and displacements caused by thermal expansion at 84 GNSS stations. (**b**,**c**) Mean precipitation in the wet (May to October) and dry (November to April) seasons from 2008 to 2012. (**d**,**e**) The monthly average precipitation and temperature in Taiwan Province, respectively.

The main source of hydrological loading is the annual rainfall in monsoon (May to June) and typhoon (July to October) seasons. A total of 70% of annual precipitation occurs during the wet season (May to October, Figure 9b), with the remaining 30% occurring during the dry season (November to April, Figure 9c). Although there is abundant annual rainfall, due to the steep terrain of Taiwan, most precipitation either evaporates directly or flows to the sea. The peak rainfall occurs from June to July (Figure 9d). Controlled by infiltration rate and capacity, the peak subsidence of GNSS lags one to two months behind precipitation, which is also consistent with the predicted value of hydrological model. At the same time, this period corresponds to the time when the highest temperatures occur in the northern hemisphere (Figure 9e). Surface water and shallow groundwater produce downward subsidence displacements on the ground, but the thermal expansion effect caused by temperature change will cause slight uplift phenomenon. This would partially counteract the observed hydrological loading displacements, and cannot be neglected in the study of water storage capacity using GNSS.

#### **5. Conclusions**

We studied the influence of non-tectonic deformation signals on the vertical coordinate time series of 84 GNSS stations in Taiwan Province. To do this, we compared the calibration effects of ATML, HYDL, and NTOL products provided by EOST, GFZ, and IMLS based on different models for Taiwan GNSS vertical time series. The differences between models were analyzed. Finally, we determined a set of model combinations suitable to correct for vertical displacements for data from GNSS stations in Taiwan Province, namely ATML (GFZ\_ECMWF IB) + HYDL (IMLS\_MERRA2) + NTOL (IMLS\_MPIOM06). We have several interesting findings from our work. Because of the special geographical location of Taiwan Province, there is an obvious monsoon alternation in Taiwan Strait, and there are many tropical cyclones from July to September every year. The tropical storms and typhoons cause abnormal changes in atmospheric and seawater movement. Because of this, non-tidal atmospheric loading and non-tidal ocean loading work together in Taiwan. Non-tidal ocean and atmospheric loading exhibit similar temporal and spatial action patterns, with better correction effect of GNSS stations in the western plain of Taiwan to that for stations in eastern Taiwan. This may be due to a central barrier of mountains as well as additional mountains in the eastern part of Taiwan.

Tropical cyclones bring abundant rainfall, which plays a vital role in groundwater recharge in Taiwan Province. However, due to the relatively small area of Taiwan, no high-resolution hydrological data model has been developed that can accurately evaluate the temporal and spatial changes of water storage capacity in Taiwan, especially due to the variations of rainfall pattern, infiltration rate, and soil saturation in temporal and spatial domains. Although existing hydrological loading data have good corrections for most GNSS stations, there are limits to the use of these data to accurately quantify hydrological changes in Taiwan Province.

We quantified the vertical displacements of different stations due to temperature changes. Because Taiwan includes middle and low latitudes, the temperature difference relative to the annual average temperature changes less than 10 ◦C, so the average annual amplitude vertical total displacement caused by thermal expansion effect is only 0.5 mm. However, maximum annual subsidence and uplift caused by hydrological loading occur in August and February respectively, with the annual phase change of displacement caused by temperature change is opposite. Thus, the counteracting effect between the thermal expansion displacement and hydrological loading displacement cannot be neglected, and the consideration of thermal expansion effect would help accurately quantify the change of water storage capacity in Taiwan Province using GNSS vertical time series.

**Author Contributions:** Conceptualization, B.L. and X.X.; methodology, B.L. and X.M.; software, B.L. and X.M.; validation, X.M.; formal analysis, B.L. and X.M.; investigation, X.M.; data curation, X.M.; writing—original draft preparation, B.L. and X.M.; writing—review and editing, J.T., W.P., and L.Z.; visualization, X.M.; supervision, B.L. and X.X.; funding acquisition, B.L. and X.X. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is funded by the National Natural Science Foundation of China (Grant No. 41904003), the Natural Science Foundation of Hunan Province, China (Grant No. 2020JJ5571), Open Fund of Engineering Laboratory of Spatial Information Technology of Highway Geological Disaster Early Warning in Hunan Province (Changsha University of Science and Technology) (Grant No. KFJ190602), Graduate research and innovation project (Changsha University of Science and Technology) (Grant No. CXCLY2022018), and Undergraduate research and innovation project (Changsha University of Science and Technology) (Grant No. 20210070).

**Data Availability Statement:** The GNSS data used in this work can be obtained from the Institute of Earth Sciences, Academia Sinica, Taiwan. The time series of the GNSS stations used in this paper are uploaded to Figshare, and readers can obtain time series data of GNSS sites classified by different monuments, as well as the vertical velocities with their uncertainties from: https: //doi.org/10.6084/m9.figshare.20448234.v3 (accessed on 7 August 2022). The environmental loads data provided by IMLS are available at http://massloading.net/#Download (accessed on 7 August 2022). The environmental loads data provided by EOST are available at http://loading.u-strasbg.fr/ displ\_maps.php (accessed on 7 August 2022). The environmental loads data provided by GFZ are available at http://esmdata.gfz-potsdam.de:8080/repository/entry/show/ (accessed on 7 August 2022). The temperature reanalysis data provided by the European Center for Medium Range Weather Forecasts are available at https://apps.ecmwf.int/datasets/data/interim-full-daily (accessed on 7 August 2022). Rainfall data provided by global precision measurement mission (GPM) comes from: https://disc.gsfc.nasa.gov/datasets?keywords=GPM&page=1 (accessed on 7 August 2022).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

