MONITORING VERTICAL LAND MOTION IN MALAYSIA USING GLOBAL POSITIONING SYSTEM (GPS) Ami Hassan Md Din12, Mohd Nadzri Md Reba2, Kamaludin Mohd Omar1, Sahrum Ses1 and Amir Sharifuddin Ab Latip3 1 Geomatic Innovation Research Group (GIG), Faculty of Geoinformation and Real Estate, Universiti Teknologi Malaysia, 81310 Johor Bahru, Johor, Malaysia. Email: [email protected], [email protected], [email protected] 2 Geoscience and Digital Earth Centre (INSTEG), Universiti Teknologi Malaysia, 81310 Johor Bahru, Johor, Malaysia. Email: [email protected] 3 Civil and Environmental Engineering Department, Universiti Teknologi PETRONAS, Perak, Malaysia. Email: [email protected] KEYWORDS: Vertical Land Motion, Global Positioning System (GPS) and Bernese ABSTRACT: The Global Positioning System (GPS) developed by the U.S. Department of Defense (DoD) for military and civilian navigation and positioning, is publicly utilized as a geodetic method of choice for investigating a wide range of geophysical phenomena. The capability of GPS in determining positions at higher accuracy has now emerged as one of the most beneficial tools for scientific research. This study intends to quantify the vertical land motion trend in Malaysia using GPS technique starting from 1999 to 2011. This study used GPS data taken from the Department of Surveying and Mapping Malaysia (DSMM) encompassing all available Continuously Operating Reference Station (CORS) data over Peninsular Malaysia, Sabah and Sarawak. A set of GPS data downloaded from the International GNSS Services (IGS) website is also used for reference. Data of 117 GPS stations were processed of which 87 stations of the Malaysia CORS networks and 30 stations of the IGS tracking network. Vertical biases of GPS measurement were reduced by applying double difference correction in processing model of Bernese v 5.0 and later transformation parameter is determined based on the minimal constraint relative to the IGS data. This processing involved epoch of 1 January 2008 and ITRF2008 reference frame. To assure the highest accuracy, the root mean square (RMS) of the single difference, the percentage of the resolved ambiguity and the RMS of daily repeatability are taken into account. The GPS derived vertical land motion trend is quantified based on the robust fit regression model. The results show that Malaysia experienced vertical land motion effects for both land uplift and subsidence of which the uplift rate ranges from 0.21 +/- 0.14 mm/yr at MUKH station to 1.44 +/- 0.13 mm/yr at PDIC station and the subsidence rate ranges from -0.04 +/- 0.04mm/yr at KUAL station to -34.41 +/- 0.16 mm/yr at AMAN station respectively. In general, land subsidence effects are more dominant. The vertical land motion was suspected to undergo local deformation as it has irregular vertical displacement patterns. However, the vertical land motion magnitude is likely to have a higher rate in the northern than in the southern part of Peninsular Malaysia. There are significant subsidence signals existed at locations such as AMAN, MUKA, KAPI, BELA, LAWS and MRDU in Sabah and Sarawak. This study concluded that the method is reliable to provide acceptable insight of the vertical land motion trend in Malaysia that may significantly benefit to other scientific works. 1. INTRODUCTION The Global Positioning System (GPS) which was developed by the U.S. Department of Defence (DoD) for military and civilian navigation and positioning, has been utilised as a geodetic method of choice for investigating a wide range of geophysical phenomena. Presently, GPS observations have been used to quantify the motion of the earth’s tectonic plates, to investigate crustal deformation around active faults and volcanoes, and to measure the displacement of the earth’s surface due to past and present changes in the mass of the world’s ice sheets (Davis et al., 2012). The capability of GPS in determining positions with high accuracy has now emerged as one of the most beneficial tools for many scientific communities. Despite the tremendous advances in GPS measurements during the last decade, a major limitation of this technique is the lack of deformation data in many areas since GPS observations are station-dependent providing only point-wise data. Currently, in Malaysia, the only active GNSS Continuously Operating Reference Stations (CORS) is the Malaysia Real Time Kinematic GNSS Network (MyRTKnet) which consists of 78 stations, with a spacing of between 30 to 100 km between one another in Peninsular Malaysia, and 30 to 200 km in Sabah and Sarawak (Mohamed, 2009). This study

intends to quantify the vertical land motion trend in Malaysia using GPS technique starting from 1999 to 2011. Vertical biases of GPS measurement were reduced by applying double difference correction in processing model of Bernese v 5.0. 2. MATERIALS AND METHODS 2.1 The GPS Data Utilised In Malaysia, the Department of Survey and Mapping Malaysia (DSMM) is responsible for providing horizontal and vertical positions for various geodetic applications. The Malaysian Active GPS System (MASS) is the first CORS network to be established in Malaysia. It was launched in early 1998. Since June 2007, the MASS CORS network has been replaced by the Malaysia Real-Time Kinematic GNSS Network (MyRTKnet). The replacement was due to the initiative of DSMM to upgrade from providing post-processed data to near real-time data (Jhonny, 2010). With the establishment of MASS and MyRTKnet CORS networks, DSMM has been continuously recording carrier phase and pseudorange measurements from all GNSS satellites. The GPS data utilised in this study was requested from the Department of Surveying and Mapping Malaysia (DSMM) encompassing all available CORS data over Peninsular Malaysia, Sabah and Sarawak. GPS data from MASS stations are available from January 1999 to December 2006 as well MyRTKnet data from December 2004 to December 2011. Figure 1 shows the distribution of CORS station used in Malaysia in this study. Besides, as for the reference station, a set of GPS data has been downloaded from the International GNSS Services (IGS) website through ftp://cddis.gsfc.nasa.gov/gps/data/daily/. 30 IGS stations around the world have been used in this study for connection of the Malaysian CORS to the International Terrestrial Reference Frame (ITRF). Figure 2 highlights the welldistributed 30 IGS stations comprising different plates. Overall, 117 GPS stations have been processed in Bernese, of which 87 stations are from MASS and MyRTKnet CORS networks and 30 stations are from IGS tracking network.

Figure 1. The distribution of CORS stations used in Malaysia 2.2 Bernese Processing Strategy for Determination of Vertical Land Motion When investigating a problem such as vertical land motion using GPS technology, the selection of processing strategy is of paramount importance. In this study, double-difference GPS processing strategy in Bernese is applied. The doubledifference technique, or normally known as relative positioning, determines the baseline vector between two stations or more if the receivers are observed simultaneously. Additionally, at least one station is held fixed (or using minimum constraint). The accuracy of the double-difference technique depends on the baseline length and satellite ephemeris accuracy which will later be used to estimate the other station coordinates (Leick, 2004).

Figure 2. Distribution of 30 IGS stations employed in this study

In Bernese, a double-difference analysis for RINEX GNSS observation data can be performed using RNX2SNX.PCF. The significant output generated from this PCF is a SINEX file that contains information about coordinates, troposphere and ambiguity parameters estimated to facilitate further processing and combination. Besides, the normal equation parameters are saved for each session to allow the estimation of station velocities. Since a large amount of GPS data that needs to be processed in this study. Hence, by applying Bernese Processing Engine (BPE) in Bernese, it has saved much processing time. In BPE, ideally, the whole process from data transfer to results presentation is automated. Bernese will run processing scripts, over multiple sessions, via the BPE. The complete processing procedure using BPE with double-difference strategy in Bernese is demonstrated in Figure 3. In data preparation, in order to obtain the highest possible precision, the IGS final satellite orbits, earth rotation parameters and absolute antenna phase centre variations and offsets were downloaded from the Bernese (http://www.bernese.unibe.ch) and IGS (http://igscb.jpl.nasa.gov/igscb/product) website. Besides, this part deals with computing a-priori coordinates, prepare local CORS and IGS observation files, create station information files, prepare ionosphere model, extrapolate station velocity file and create ocean loading data. Ocean tide loading parameters can be downloaded from the free service website provided at http://www.oso.chalmers.se/~loading/. Next stage is data transfer. This part of the script deals with converting RINEX files to Bernese format, synchronising the receiver clocks to GPS time, and producing a summary of available data and their performance. Two Bernese programs are involved in this step: RXOBV3 and RNXGRA. Then, the information such as pole, orbit and clock need to be prepared. Three Bernese programs are used in this part. They are POLUPD to generate the pole information file in Bernese Format (Orbit and Earth orientation files are converted from foreign to Bernese format), PRETAB to generate tabular orbit files from IGS precise orbit (conversion of precise orbit information from SP3 into the inertial frame / tabular format), and ORBGEN to generate Bernese standard orbit files (the tabular orbit file is used as input for the final orbit generation step in ORBGEN. The result is a standard orbit file and a summary providing the quality of the orbit fit).

Figure 3. GPS double-difference processing flow in Bernese using BPE

The screen code data part deals with computing the corrections for synchronising receiver clock with respect to GPS time (i.e. receiver clock corrections) by using program CODSPP. Only code observations are used for this step. The receiver clock corrections computed by CODSPP are stored in CODE+PHASE observation file. In From Baselines stage, Program SNGDIF forms phase single-difference baselines from zero-difference observation files. The adopted strategy for the selection process is OBS-MAX. This strategy guarantees that a baseline is formed by taking into account the maximum number of observations between stations to form a network for the double-difference processing. The procedure is; however, time consuming because the zero-difference observation files have to be analysed several times. Next, Program MAUPRP detects and resolves cycle slips, removes outliers, and adds multiple ambiguities for the phase observation files in Screen Phase Data stage. It works with station observation files (zero-difference mode) and with baseline observation files (single-difference mode). If the size of a cycle slip cannot reliably be determined, a new ambiguity is set up. In Baseline Processing, this step deals with computing ambiguity-FLOAT network solution and resolve phase ambiguity. Two Bernese programs are involved: GPSQIF and ADDNEQ2. Program GPSQIF, the Quasi Ionosphere Free (QIF) strategy, is used to resolve ambiguities to integer number. In this step, the processing will be made for each baseline separately. While for program ADDNEQ2, a network solution with real valued ambiguities is computed based on normal equations and stored in GPSEST after residual screening. “Float” coordinates are saved for further use in the ambiguity resolution step. The final part is to compute the final ambiguity fixed solution after all the observation files are cleaned and most of the ambiguities are resolved to their integer values. The final results use the programs GPSEST and ADDNEQ2. Program GPSEST computes an ambiguity fixed solution and normal equation information is stored. Estimated parameters include station coordinates, zenith path delays and horizontal troposphere gradients. On the other hand, ADDNEQ2 computes final coordinate solution based on the normal equations from the previous GPSEST. The datum definition is realized by three no-net-translation, three no-net rotations and scale conditions imposed on a set of fiducial sites (ITRF2008 reference station coordinates). The transformation parameter computation is based on the minimal constraint of the IGS station. The advantages of a minimum constraint solution are the minimal errors in the coordinate of a reference site and the minimal distortion of network geometry (Dach et al., 2007). Table 1 shows the summarised processing parameters and models that have been applied to obtain high quality vertical estimates in Bernese software version 5.0 using BPE. Table 1. Processing parameters and models for GPS data processing Processing Parameters Processing Strategy Daily Input data OBS-MAX Network design 3° Elevation cut off angle 30 s Sampling rate IGS final orbits (SP3) and EOP (Earth Orientation Parameters). Orbits / EOP Minimally constrained to the ITRF2008 reference frame. Reference Frame PHAS COD.I08 Antenna phase centre FES2004 Ocean loading model Double difference ionospheric free (IF) linear combination (L3). Ionosphere Fixed, resolved using QIF strategy strategy with baseline < 2000km Ambiguities solution Ionosphere model for Global ionosphere model from CODE ambiguity fixing A-priori Saastamoinen model (hydrostatic part) with dry Niell A-priori model mapping function. Wet-Niell mapping function (2h interval) Mapping Function Horizontal gradient parameters: tilting (24 h interval). Gradient Estimation

2.3 Time Series Analysis of Vertical Land Motion using Robust Fit Technique The time series of the vertical land motion trend in this study was quantified using robust fit regression analysis. Robust fit analysis is a standard statistical technique that concurrently deals with solution determination and outliers detection. In this robust fit approach, a linear trend is fitted to the annual sea level time series of each station in an Iteratively Reweighted Least Squares (IRLS) technique (Holland and Welsch, 1977). Depending on the deviations from the trend line, weights of measurements are adjusted accordingly. The trend line is then re-fitted. The process is repeated until the solution converges. The weights of the observations (wi) are readjusted by the adopted bi-square weight function, whose relationship with normalised residuals, (ui) can be written as (Holland and Welsch, 1977):

wi = {

(

( ) )

(1)

Where, ui =

√

ri : Residuals, hi : Leverage, S : Mean absolute deviation divided by a factor 0.6745 to make it an unbiased estimator of standard deviation K : A tuning constant whose default value of 4.685 provides for 95% asymptotic efficiency as the ordinary least squares assuming Gaussian distribution Observations that are assigned zero weights in any iteration are declared as outliers and eliminated from further computation (Holland and Welsch, 1977). 3. RESULTS AND DISCUSSION 3.1 GPS Data Quality Control For GPS data quality control particularly in the program GPSQIF, there are two outputs that need to be taken into account. The first is the value of root mean square (RMS) of the single difference and the second one is the percentage of the resolved ambiguity. RMS of a single-difference baseline for good GPS data is smaller than 2 mm (Dach et al., 2007). The resolved ambiguity percentage for good GPS data should be higher than 75 %. Table 2 shows the ambiguity resolution results for Day of Year (DOY) 30, 2010 from program GPSQIF. The ambiguity resolution for this sample has shown good results as its average is 84.9% and 1.7 mm for RMS of single difference. Table 2. Good ambiguity resolution summary (DOY 30, 2010 data)

In addition, it is necessary to study the RMS of repeatability of ITRF2008 mapped daily solutions in order to understand the ability and sensitivity of GPS data in vertical deformation detection. In GPS measurements, the precision of a point’s coordinate is defined as the formal standard error associated with the individual determination, and is normally calculated as the Root Mean Square (RMS) (Hofton, 1995; Williams, 1995). Figure 4 shows the RMS of daily repeatability with respect to (w.r.t) monthly averaged solutions for GETI, KUAL, MIRI, MTAW, SAND and USMP stations (sample station). The RMS of daily repeatability throughout the whole GPS data set typically range from 0.3 to 1.0 mm for horizontal component (north and east) and 0.8 to 1.8 mm for vertical component (up). The RMS of horizontal component normally is three times better than the RMS of vertical component. From Figure 4, the RMS results are quite interesting, where the RMS of daily repeatability for the year 2005 and onwards depict better results than the years before. The RMS value for the year 2005 and onwards typically range from 0.3 to 0.4 mm for horizontal component and 0.8 to 1.2 mm for vertical component. This is possibly due to the GPS stations’ good performance after new hardware installation, as it is during the transition period from MASS stations to the MyRTKnet stations. In summary, the RMS of daily repeatability for all GPS stations are considered excellent and also proof that GPS data processing via Bernese software performed excellently well in this study in deriving horizontal and vertical time series coordinates.

Figure 4. Daily repeatability w.r.t monthly averaged solutions for GETI, KUAL, MIRI, MTAW, SAND and USMP stations 3.2 GPS-derived Vertical Land Motion Rate The GPS-derived vertical land motion was produced by the high precision Bernese software using double differenced strategy. In order to quantify the vertical displacement, the time series of station coordinates are referred to ITRF2008 at epoch 01 January 2008. The GPS-derived vertical land motion rate is quantified using robust fit regression technique in MATLAB. The GPS data time frame is analysed from 1999 up to 2011, but most newly established MyRTKnet stations only commenced after year 2004. Hence, the time frame for GPS data analysis is dissimilar between one another for most GPS stations in this study.

Since the aim of this paper is to quantify the rate of vertical land motion over the Malaysian region, only the vertical component from the GPS results are analysed in this section. Figure 5 demonstrates the plot of vertical displacement time series analysis in daily solutions using robust fit regression at GETI, KUAL, MIRI, MTAW, SAND and USMP stations (sample station). The red dotted line represents the original (observed) daily GPS data for the vertical component, the blue line illustrates the robust fit simulation line after applying iteratively re-weighted least square (IRLS) and the black dotted line represents the linear trend using robust fit regression. As expected, from Figure 5, the daily consistency for vertical component shows better results from the year 2005 and onwards. These vertical displacement consistency results are tallied exactly with the daily precision results that have been discussed previously in Section 3.1. For some stations, for instance GETI, the annual oscillation is clearly seen in the daily vertical displacement time series. This displacement cycle is caused by seasonal changes in atmospheric pressure.

Figure 5. Vertical displacement time series in daily solutions for GETI, KUAL, MIRI, MTAW, SAND and USMP stations Table 3 presents the summary of GPS-derived vertical land motion rate and their uncertainties (standard errors) in mm/yr over the Malaysian region. From Table 3, the Malaysian region has vertical land motion effects for both land uplift and subsidence. The uplift rate ranges from 0.21 +/- 0.14 mm/yr at MUKH to 1.44 +/- 0.13 mm/yr at PDIC. On the other hand, the rate of subsidence over the Malaysian region ranges at rate of -0.04 +/- 0.04mm/yr at KUAL to 34.41 +/- 0.16 mm/yr at AMAN. Interestingly, there are three GPS stations that have subsidence rates more than 1 cm/yr. They are LAWS at a rate of -13.83 +/- 0.24 mm/yr, MUKA at a rate of -17.47 +/- 0.25 mm/yr and AMAN at a rate of -34.41 +/- 0.16 mm/yr. This may be due to the fact that those affected areas are located at valley areas. Hence,

these stations tend to subside easily. However, further investigation needs to be carried out at those affected areas in order to know the true source of these significant subsidence rates. Table 3. The GPS-derived vertical land motion rates and their uncertainties (standard errors) in mm/yr over the Malaysian region derived from Bernese software No 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45

Stn Name AMAN ARAU AYER BABH BAHA BANT BEAU BEHR BELA BELU BENT BIN1 BINT CAME CENE DATU GAJA GETI GMUS GRIK IPOH JHJY JRNT JUML KAPI KENI KINA KLAW KRAI KROM KTPK KUAL KUAN KUCH KUDA KUKP LAB1 LABU LASA LAWS LGKW LIPI MERS MERU MIRI

Latitude 1°13'18.70" 6°27'00.57" 5°45'00.88" 5°08'47.97" 2°48'23.42" 2°49'33.44" 5°21'21.09" 3°45'55.33" 3°02'43.14" 6°25'02.84" 3°31'36.91" 3°14'25.14" 3°15'41.42" 4°25'28.99" 4°07'23.14" 5°01'38.35" 2°07'20.24" 6°13'34.29" 4°51'46.71" 5°26'20.44" 4°35'18.50" 1°32'12.52" 3°55'27.41" 2°12'42.31" 2°00'26.67" 5°22'03.19" 5°54'16.57" 2°58'53.43" 5°30'07.17" 2°45'47.02" 3°10'15.40" 5°19'08.00" 3°50'03.76" 1°37'56.64" 6°53'53.31" 1°19'59.79" 5°16'57.51" 5°16'57.61" 4°55'25.81" 4°52'19.86" 6°19'42.60" 4°10'33.62" 2°27'12.48" 3°08'17.65" 4°22'19.56"

Longitude 111°27'27.24" 100°16'47.06" 101°51'36.53" 100°29'37.18" 102°22'40.37" 101°32'14.47" 115°43'49.48" 101°31'01.96" 113°54'57.44" 116°26'43.57" 101°54'25.93" 113°05'39.61" 113°04'02.00" 101°23'07.07" 103°14'23.17" 118°17'49.47" 103°25'21.76" 102°06'19.67" 101°57'49.66" 101°07'48.99" 101°07'34.24" 103°47'47.52" 102°23'03.89" 102°15'21.96" 112°55'40.77" 116°10'51.35" 116°02'21.48" 102°03'49.20" 102°13'10.86" 103°29'50.27" 101°43'03.39" 103°08'20.93" 103°21'01.27" 110°11'42.33" 116°50'54.50" 103°27'12.36" 115°14'41.19" 115°14'41.24" 101°04'04.95" 115°24'45.47" 099°51'04.54" 102°06'01.94" 103°49'43.51" 101°24'26.85" 114°00'06.28"

Data Span 092007 - 122011 011999 - 122011 072007 - 122011 122004 - 122011 062007 - 122011 122004 - 122011 122007 - 122011 122002 - 122011 112008 - 122011 112008 - 122011 062007 - 122011 092007 - 122011 021999 - 102006 062007 - 122011 062007 - 122011 112008 - 122011 072007 - 122011 011999 - 122011 122004 - 122011 122004 - 122011 041999 - 082006 122004 - 122011 072007 - 122011 122004 -122011 122007 - 122011 122007 - 122011 011999 - 122006 122004 - 122011 072007 - 122011 062007 - 122011 041999 - 122006 011999 - 122011 031999 - 122006 011999 - 112006 122007 - 122011 122004 - 122011 032006 - 122011 011999 - 012006 072007 - 122011 112008 - 122011 122004 - 122011 072007 - 122011 122004 - 122011 122004 - 122011 031999 - 122011

VLM Rate (mm/yr) -34.41 +/- 0.16 -1.17 +/- 0.05 -1.34 +/- 0.14 -3.17 +/- 0.07 0.38 +/- 0.12 -5.98 +/- 0.17 -1.64 +/- 0.16 -0.95 +/- 0.06 -5.07 +/- 0.51 -1.28 +/- 0.24 0.24 +/- 0.13 -2.62 +/- 0.16 -0.23 +/- 0.11 -0.28 +/- 0.14 -0.81 +/- 0.13 -0.23 +/- 0.26 -1.96 +/- 0.13 -1.98 +/- 0.04 -2.50 +/- 0.08 -2.65 +/- 0.08 -1.62 +/- 0.10 -2.12 +/- 0.07 0.91 +/- 0.12 -0.98 +/- 0.06 -6.99 +/- 0.15 -1.00 +/- 0.17 -0.22 +/- 0.10 -1.57 +/- 0.07 -6.09 +/- 0.15 -2.04 +/- 0.12 -0.26 +/- 0.08 -0.04 +/- 0.04 -0.68 +/- 0.08 -1.20 +/- 0.14 -1.05 +/- 0.14 -5.10 +/- 0.06 -0.30 +/- 0.14 -2.48 +/- 0.10 0.48 +/- 0.12 -13.83 +/- 0.24 -3.53 +/- 0.07 -6.00 +/- 0.12 -2.86 +/- 0.07 -8.92 +/- 0.08 0.51 +/- 0.04

46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87

MRDI MRDU MTAW MUAD MUKA MUKH NIAH PASP PDIC PEKN PRTS PUPK PUSI RANA SAND SARA SBKB SEG1 SEGA SEMA SEMP SETI SGPT SIB1 SIK1 SPGR SRIJ TEBE TENM TERI TGPG TGRH TLKI TLOH TMBN TOKA UMAS UMSS UPMS USMP UTMJ UUMK

4°11'01.46" 6°31'26.01" 4°15'45.98" 3°04'18.45" 2°52'21.66" 4°37'03.49" 3°51'43.94" 5°50'16.63" 2°31'34.23" 3°29'33.35" 1°58'53.07" 4°12'25.17" 4°28'50.52" 5°58'25.81" 5°50'32.65" 1°45'50.78" 3°48'45.99" 2°29'10.68" 2°29'10.67" 1°48'09.69" 4°28'16.38" 5°31'56.98" 5°38'36.88" 2°04'18.91" 5°48'35.64" 1°48'38.14" 3°39'35.77" 1°01'01.06" 5°10'05.45" 5°08'49.11" 1°22'02.68" 2°04'46.76" 3°59'28.80" 3°26'58.02" 5°43'25.89" 6°01'46.59" 1°28'05.91" 6°02'21.24" 2°59'36.22" 5°21'28.03" 1°33'56.93" 6°27'43.85"

114°20'01.47" 116°46'09.73" 117°52'53.94" 103°04'27.98" 112°01'11.68" 103°12'34.02" 113°42'51.84" 102°21'27.62" 101°48'37.92" 103°23'22.89" 102°52'23.03" 100°33'33.28" 101°01'06.34" 116°40'29.49" 118°07'14.11" 111°20'11.69" 100°48'59.06" 102°43'55.27" 102°43'55.35" 109°45'48.51" 118°37'05.36" 102°43'57.30" 100°29'18.15" 111°40'23.76" 100°43'44.01" 103°19'15.53" 102°54'30.16" 110°21'18.49" 115°57'37.76" 102°57'54.43" 104°06'29.74" 103°56'48.99" 101°03'13.82" 102°25'09.72" 116°24'10.61" 100°24'12.86" 110°25'28.92" 116°06'43.31" 101°43'24.64" 100°18'14.54" 103°38'22.43" 100°30'22.81"

112008 - 122011 112008 - 122011 061999 - 122011 062007 - 122011 022009 - 122011 072007 - 122011 112008 - 122011 072007 - 122011 062007 - 122011 122004 - 122011 062007 - 122011 122004 - 122011 072007 - 122011 112008 - 122011 051999 - 122011 112008 - 122011 062007 - 122011 072007 - 122011 022001 - 112006 122007 - 122011 112008 - 122011 062007 - 122011 122004 - 122011 122007 - 122011 072007 - 122011 072007 - 122011 062007 - 122011 112008 - 122011 112008 - 122011 072007 - 122011 122004 - 122011 062007 - 122011 072007 - 122011 122004 - 122011 122008 - 122011 042008 - 122011 122004 - 122011 032005 - 122011 122004 - 122011 121999 - 122011 011999 - 122006 122004 - 122011

-0.82 +/-9.90 +/-0.09 +/-0.66 +/-17.47 +/0.21 +/-4.57 +/-1.11 +/1.44 +/-3.30 +/-1.35 +/-0.98 +/-0.40 +/-0.50 +/-1.58 +/-0.92 +/-0.39 +/-2.18 +/-1.39 +/-0.31 +/-0.47 +/-1.36 +/-2.52 +/-0.78 +/-0.67 +/-0.14 +/-7.03 +/-2.43 +/-0.12 +/-2.20 +/-3.43 +/-1.15 +/-0.68 +/-1.05 +/-0.62 +/-5.80 +/-0.53 +/-3.69 +/-0.89 +/-0.90 +/-0.62 +/-4.38 +/-

0.24 0.28 0.04 0.11 0.25 0.14 0.21 0.14 0.13 0.06 0.13 0.06 0.12 0.25 0.04 0.23 0.11 0.14 0.14 0.17 0.23 0.13 0.07 0.15 0.14 0.12 0.13 0.25 0.31 0.14 0.06 0.12 0.13 0.06 0.26 0.17 0.06 0.08 0.06 0.04 0.09 0.07

Using colours on a map is a practical method to visualise vertical displacements over the Malaysian region. Figure 6 demonstrates the GPS-derived vertical displacement trend map over Peninsular Malaysia, Sabah and Sarawak. From the figure 6, land subsidence effects are dominant over the Malaysian region. The vertical land motion seems to undergo local deformation as it has irregular patterns in vertical displacement. However, the vertical land motion magnitude seems to have a higher rate in the northern part of Peninsular Malaysia comparatively to the southern part of Peninsular Malaysia. While in Sabah and Sarawak, there are significant subsidence signals at a few (station) locations such as AMAN, MUKA, KAPI, BELA, LAWS and MRDU. The vertical land motion is less constrained than the horizontal motion due to the uncertainties of GPS measurements.

Figure 6. Vertical land motion trend colour map derived from GPS data over Peninsular Malaysia and, Sabah and Sarawak. Units are in mm/yr. 4. CONCLUSIONS The results show that Malaysia experienced vertical land motion effects for both land uplift and subsidence of which the uplift rate ranges from 0.21 +/- 0.14 mm/yr at MUKH station to 1.44 +/- 0.13 mm/yr at PDIC station and the subsidence rate ranges from -0.04 +/- 0.04mm/yr at KUAL station to -34.41 +/- 0.16 mm/yr at AMAN station respectively. In general, land subsidence effects are more dominant. This study concluded that the method is reliable to provide acceptable insight of the vertical land motion trend in Malaysia that may significantly benefit to other scientific works. ACKNOWLEDGEMENTS The authors would like to thank to Department of Survey and Mapping Malaysia (DSMM) for providing GPS data. We are grateful to the Ministry of Science, Technology and Innovation (MOSTI) for funding this project under the eScience Fund, Vote Number 04-01-06-SF1092. REFERENCES Dach, R., Hugentobler, U., Fridez, P. and Meindl, M. (2007). Bernese GPS Software Version 5.0. Astronomical Institute, University of Bern, http://www.bernese.unibe.ch/docs/DOCU50.pdf. Davis, J. L., Prescott, W. H., Svarc, J. L. and Wendt, K. J. (2012). Assessment of Global Positioning System Measurements for Studies of Crustal Deformation. Journal of Geophysical Research: Solid Earth (1978–2012), Vol 94 Issue B10 pages 13635–13650. DOI: 10.1029/JB094iB10p13635. Hofton, M. A. (1995). Anelastic Deformation in Iceland Studied using GPS: with Special ReferenceTtopost-tectonic Motion following the 1975-1985 Krafla Rifting Episode, and Isostatic Rebound. Doctor Philosophy, Durham University, UK. Holland, P. W. and Welsch, R. E. (1977). Robust Regression using Iteratively Reweighted Least-squares. Communications in Statistics—Theory and Methods 6 (9), 813–827.

Jhonny (2010). Post-seismic Earthquake Deformation Monitoring in Peninsular Malaysia using Global Positioning System. MSc Thesis, Universiti Teknologi Malaysia, Skudai. Leick, A. (2004). GPS Satellite Surveying. (3rd ed.). New Jersey:John Wiley and Sons Inc. Mohamed, A. (2009). JUPEM GNSS Infrastructure. Workshop on Surveying with a Single GPS Receiver in MyRTKnet Enviroment, 29-31 July. Kuching – Sarawak. Williams, S. D. P. (1995). Current Motion on Faults of the San Andreas System in Central California Inferred from Recent GPS and Terrestrial Survey measurements. Doctor Philosophy, Durham University, UK.