c ESA/ESTEC, 2012 20 years of Progress in Radar Altimetry, 24-26 September 2012, Venice, Italy, vol. ESA SP-710.

SEA LEVEL RISE AND SUBSIDENCE IN THE DELTA AREAS OF THE GULF OF THAILAND Marc Naeije1 , Wim Simons1 , Itthi Trisirisatayawong2 , Chalermchon Satirapod2 , and Sommart Niemnil3 1

Aerospace Engineering, TUDelft, Kluyverweg 1, 2629 HS Delft, The Netherlands, Email:[email protected] 2 Dept. Survey Engineering, Chulalongkorn University, Patumone Bangkok 10330, Thailand 3 Royal Thai Naval Academy, 204 Sukumvith Rd. Paknam, Samutprakan Bangkok 10270, Thailand

ABSTRACT In the Thailand-EC GEO2TECDI projects we investigate the vertical land motion in Thailand and sea level change in the Gulf of Thailand. Focus is on Bangkok which 1) is situated in a river delta and average height is close to sea level, 2) is subsiding due to ground water extraction, 3) is experiencing post-seismic motion due to nearby megathrust Earthquakes, and 4) suffers from rising sea levels due to global climate change. This poses a serious threat to Thai society and economy. Before mitigation methods can be devised we aim at charting, qualifying and quantifying all contributing effects by the use of GNSS, InSAR and altimetry, and combining results with in situ observations like tide gauges and with geophysical modelling. Adding GPS based vertical land motion to the tide gauge sea level registration reveals the absolute sea level change at a number of tide gauge stations, which is nicely confirmed by altimetry. We find an average absolute rise of ≈3.5mm/yr ± 0.7, but near estuaries of the Chao Praya River (Bangkok), the Kah Bpow River (Koh Kong) and the Mekong delta (Ho Chi Min City), this mounts to 4 to 5mm/yr, faster than the global average. This is reinforced when accounting for the tectonic subsidence that resulted from the 2004 9.1Mw Sumatra/Andaman earthquake; from 2005 onwards we find downfall in the order of 10mm/yr resulting in relative sea level rise of >10mm/yr over the last 7 years. RADARSAT InSAR analyses confirm this and show subsidence rates up to 25mm/yr at many places along coastal Bangkok. Key words: regional sea level change, vertical land motion, subsidence, GNSS, InSAR, altimetry, tide gauges.

1.

INTRODUCTION

The EU funded GEO2TECDI projects aim to enhance skills and Thailand-EU joint research capabilities of Thai researchers in exploiting space geodetic techniques to detect, monitor and model land surface (vertical) motion due to tectonic (earthquake) processes, (human induced) land subsidence and (climate change related) sea level

changes. GEO2TECDI results already indicated that the significant magnitudes of these 3 phenomena in Thailand are strongly correlated and in the long term could have a considerable impact on the social-economic development of coastal and low-lying areas (like the capital city of Bangkok). Therefore we strive to extend, refine and improve these findings with a clear focus on the problem of (relative) sea level change threatening Thai coastal regions. This is combined with a distinct outreach to the possibly affected Thai public and the (local) authorities involved in coastal management, to timely mitigate climate change effects, to adapt or at least to be prepared.

This paper addresses the absolute and relative sea level change in the Gulf of Thailand combining the 20 year historic satellite altimeter data record with vertical land motion from the Global Navigation Satellite System GNSS and in situ data from local tide gauges. Studying sea level change rate from tidal data is complicated by the fact that tide gauges are attached to the land and thus their measurements are affected by vertical land motion caused by natural or anthropogenic processes unrelated to variations in the absolute sea level. In the traditional approach of many global studies, subsets of tide gauges have been selected away from tectonic active areas. Introducing the land effect of ice melt since the last ice age (global isostatic adjustment) a record of the last 60 years gives an estimate of absolute rates of sea level change at around 1.8 mm/yr for the 20th century. This by no means is a good estimate for regional or local sea level change, so also not for the Gulf of Thailand. Though GIA may be neglected here, land deformation due to tectonic processes not. This might appear strange because Thailand itself was never hit by a major earthquake, but it still suffers the consequences from the large earthquakes in the Indonesian Archipelago, especially the Sumatra/Andaman 9.1Mw earthquake at the end of 2004. Even in Bangkok, more than 1000km away from the epicentre we notice the coseismic motion, and more importantly the post-seismic motion. This is exactly what we want to address in this paper. In Section 2 we explain the data and methodologies applied, in Section 3 we present and discuss our results, and we end the paper with some conclusions in Section 4. Page 1 of 6

c ESA/ESTEC, 2012 20 years of Progress in Radar Altimetry, 24-26 September 2012, Venice, Italy, vol. ESA SP-710. As the quality of the radial accuracy of the orbits of the TOPEX-class satellites (TOPEX, POSEIDON and JASON) outperforms the accuracy of the orbits of the ERSclass satellites (ERS and ENVISAT), basically because of the orbital height and influence of atmospheric drag, we performed a so-called dual-crossover minimization analysis in which the orbit of the TOPEX-class satellites is held fixed and those of the ERS-class satellites were adjusted simultaneously. Then the ERS-class and TOPEXclass data were merged and gridded to daily sea level anomaly grids. Our study area is 5◦ –14◦ North, and 99◦ – 105◦ East. We choose a rectangular mesh with block size 0.25◦ and an e-folding of 0.75◦ and cut-off at 1.5◦ (space) and of 9 days and cut-off at 18 days (time). Subsequently the daily solutions were collected per month and a monthly average was calculated. This method aims to make the final monthly altimeter solution comparable with the monthly tide gauge solution. It appeared that this approach improves the root mean square of the difference between monthly-solutions of altimeter and tide gauge (RMSE), and also the correlation. Taking the monthly mesh solutions we can determine sea level trend on a mesh point to mesh point basis. For each mesh point the time series of sea level anomaly has been subjected to a robust regression analysis using iterative reweighed least squares IRLS [4], which is also used to estimate sea level trends from the tide gauge data. Figure 1. Geographical distribution of the data sets used: Altimetry, tide gauge and GPS. 2.

DATA AND METHODOLOGIES

Figure 1 introduces the location of the data used in our study. The blue and grey lines represent the typical sampling from respectively the satellite altimeter reference 10-day repeat track by TOPEX (TOPEX/Poseidon and the two Jasons) and the higher spatial density reference 35-day repeat track by Envisat (both ERS’es and Envisat). The circles represent tide gauge locations in various flavours: colours indicating long (dark blue) and short records (light blue), and those not tied to the bedrock (white) or data unavailable (orange). The triangles represent campaign GPS locations (red), and continuous GPS locations (black). For the delta area around Bangkok we also use RADARSAT SAR images in our InSAR processing scheme to extract land motion. Altimetry To obtain altimeter data we use the Radar Altimeter Database System RADS. This database provides a continuous and consistent observing system and has the complete backlog of all satellite altimeter observations [2]. A satellite altimeter functions like a spaceborne tide gauge obtaining absolute sea level unaffected by land motion. RADS takes care of all the necessary corrections and models. Still detailed inspection is needed because we are analysing data close to the coast. In [4] a summary is given of the chosen corrections. Altimeter data has been used from ERS-1, ERS-2, Envisat, TOPEX, Poseidon, Jason-1, and Jason-2, starting January 1993 and ending January 2010.

GNSS To determine vertical land motion in Thailand we made use of the dual frequency GPS data from several GPS campaign sites and continuous GPS (CGPS) stations (see Figure 1). We have adopted the absolute GPS positioning technique referred to as precise point positioning (PPP). With this technique a daily repeatability is to be expected of a few mm horizontally and about a cm vertically. The data have been uniformly processed with the GIPSY-OASIS II software developed by the JPL. Details on this processing can be found in [3]. We have 2 stations in the vicinity of the used tide gauges: CHON 26km away from Ko Sichang and 56km from Sattahip, and BAHN 27km away from Ko Mattaphon. We also included the tide gauge at Ko Lak in our analysis because of the good fit with altimetry. For Geting (Malaysian tide gauge) we take the GPS result from [5]. The 1994–2004 station velocities were extracted by robust linear regression, implying a steady-state motion, which clearly not holds after the Sumatra 2004 9.1Mw earthquake. The Sumatra-Andaman earthquake completely changed the crustal motion pattern in Thailand, in both horizontal and vertical direction. GPS data from campaigns and fixed (continuous) sites in Indonesia, Malaysia and Thailand show that the Sundaland block, on which all tidal stations in Thailand are situated, was considerably deformed. Data from yearly GPS campaigns (unlike continuous GPS) after the earthquake are less suitable for the accurate detection of vertical co- seismic displacements, since they typically have noisier position estimates, and also are affected by a strong post-seismic signal. This makes the use of GPS data from post-2004 campaigns at BANH and CHON less reliable and measurements from continuous GPS (CGPS) more appropriate. Page 2 of 6

c ESA/ESTEC, 2012 20 years of Progress in Radar Altimetry, 24-26 September 2012, Venice, Italy, vol. ESA SP-710. MEDIUM TERM TREND 1993-2004

LONG TERM TREND 1993-2009

TREND ERROR 1993-2009

Figure 2. Geographical distribution of sea level trend and trend error from multi-satellite altimetry Geophysical modelling Since post-seismic motions are not steady rates, interpretation of tide gauge data after the earthquake is not as straightforward as in the pre-earthquake period and longer time series are required before more reliable post-seismic velocities can be obtained. In the long run the cumulative effect of postseismic tectonic plate subsidence will make it imperative for sea level change study to know these subsidence rates and the changes in these rates (either accelerations or decelerations). In order to quantify the role of asthenospheric relaxation for both the horizontal and vertical post-seismic deformations observed after the 2004 Sumatra earthquake, we have adopted a 3D finite element method [1], and compare the results obtained from our preferred set of rheological parameters with the subsidence we observe in GPS and tide gauges. We compute the surface deformation induced by the visco-elastic relaxation in the low viscosity regions of our model. The observed GPS horizontal and vertical post-seismic displacements are interpreted as a consequence of sliding on the subduction plane and relaxation in the asthenosphere and low viscosity wedge [1]. InSAR InSAR and GPS observations are two complementary geodetic techniques. Whereas GPS gives absolute 3D information on surface deformation at single locations, inSAR gives a high-resolution image of the surface displacement but inherently suffers from a systematic error (phase plane). Basically the combination of the two techniques would enable the elimination of this error, though here we are not interested in absolute displacements but in displacement rates. Another issue that needs attention is the effect of the activity of the atmosphere which generates noise errors in the observations. For this a vectorized filtering is used to reduce the noise of the modulated (wrapped) interferograms. Our

interferometry analysis is based on so-called persistent scatterers (PS), and we make use of both the Stanford method for persistent scatterers (StaMPS) and the Delft PS-InSAR processing package (DePSI). In the Bangkok area we dispose of 13 scenes ALOS-PALSAR (path 487, frame 260) from 2007 to 2010, 23 scenes RADARSAT1 (path 2, row 62 and 64) from 2005 to 2012, and 4 scenes RADARSAT-2 for 2010 and 2011. Clearly, the RADARSAT-1 data is most suited for credible land subsidence estimates, though sketches the situation only after 2005. We have processed these data with StaMPS.

3.

RESULTS AND DISCUSSION

In Figure 2 we show the geographical distribution of sea level trend from 1993–2004 (left) and 1993–2009 (centre) multi-satellite altimetry. A simultaneous fit of linear trend, annual cycle, semi-annual cycle (very small), and a bias, subjected to IRLS was applied to monthly-averaged altimeter meshes. The right panel gives the error associated to the long term trend estimate. In a direct comparison with tide gauge and GPS data we have to take into account the discontinuity in land motion due to the 2004 Sumatra earthquake. Though the sea level trend over the medium length period appears significantly affected by the lunar nodal tide effect (18.6 year cycle), the effect in tide gauge data is the same and cancels out when the two data types are differenced. We find an average absolute rise of 3.5 mm/yr ± 0.7. However, near the estuaries of large rivers like the Chao Praya River (Bangkok), the Kah Bpow River (Koh Kong) and the Mekong delta (Ho Chi Min City), this mounts to 4 to 5mm/yr, all faster than the global average. Page 3 of 6

c ESA/ESTEC, 2012 20 years of Progress in Radar Altimetry, 24-26 September 2012, Venice, Italy, vol. ESA SP-710.

Figure 3. Vertical land motions at CHON and BANH from seasonally corrected precise GPS height coordinates The vertical land motions at CHON and BANH based on GPS campaign data have been plotted in Figure 3. Heights have been corrected for seasonal variations estimated from nearby CGPS data. The dotted vertical line marks the time of the Sumatra earthquake. We see a combination of interseismicity and GIA uplift at both locations before 2005 and post-seismicity after the earthquake at the end of 2004. The difference in vertical motion regimes before and after the quake are striking (opposite and larger rates after 2004) and the post-seismicity clearly adds to the flooding threat of the low-lying parts of Thailand. All GPS campaigns were conducted by the Royal Thai Survey Department and analysed in cooperation with TUDelft. Simons et al. [3] provide details on the computational process and the 1σ error bars. Assuming that both tide gauges and altimetry capture the same absolute sea level signals, these cancel when the tide gauge trend is subtracted from the corresponding altimetry trend. The remainder then, besides measurement errors, represents the vertical land motion. Definitely, first a good match between altimetry and tide gauge time series has to be achieved, which is monitored by corre-

Tidal Station Laem Singh Sattahip Ko Sichang Hua Hin Ko Lak Ko Mattaphon Ko Prab Songkhla Geting

lation and RMSE. The tide gauge monthly-average time series and those from altimetry, both from 1993–2004 (prior to 2004 Sumatra Earthquake), have been compared and differenced for 9 stations. The intra-annual, seasonal and inter-annual signals compare very well: RMSEs from 5cm to 10cm and correlations from 0.85 to 0.96, depending on the distance of the closest TOPEX-class track passage to the tide gauge (turning point at ≈100km). Table 1 provides a summary of the altimeter and tide gauge based sea level trends over the same time period. In the 5th column the tide solution has been subtracted from the altimeter to derive the local vertical land motion (VLM). We compare this with VLM based on GPS in the rightmost column. The stations with the best fit between altimetry and tide gauge allow for an independent measure of VLM which is confirmed by GPS all within the indicated error margins. This also shows the need of collocated GPS at the exact locations of the tide gauge. This was brought to the attention of the Thai authorities. Ko Sichang is considered an outlier for a number of reasons: shorter record, low correlation and high RMSE, and a chance being tied to the sediment rather than the bedrock.

Altimetry 1993-2009

Altimetry 1993-2004

Tide Gauge 1993-2004

VLM (Alt−TG) 1993-2004

3.6 ± 0.7 4.8 ± 0.7 5.8 ± 0.8 4.8 ± 0.8 3.6 ± 0.7 3.2 ± 0.7 3.3 ± 0.6 3.3 ± 0.6 3.9 ± 0.6

6.6 ± 1.2 6.9 ± 1.1 7.8 ± 1.3 6.4 ± 1.2 5.5 ± 1.0 5.8 ± 0.9 6.0 ± 1.0 4.6 ± 1.0 6.1 ± 0.8

10.5 ± 1.1 3.2 ± 1.6 7.6 ± 2.0 4.8 ± 2.2 1.9 ± 1.1 5.8 ± 1.9 9.9 ± 1.9 13.0 ± 2.4 2.6 ± 1.2

−3.9 ± 1.6 3.7 ± 1.9 0.2 ± 2.4 1.6 ± 2.5 3.6 ± 1.5 0.0 ± 2.1 −3.9 ± 2.1 −8.4 ± 2.6 3.5 ± 1.4

VLM (GNSS) 1994-2004 N/A 3.8 ± 1.3a 3.8 ± 1.3a N/A N/A 2.2 ± 0.8b N/A N/A 5.1 ± 1.0c

Table 1. Altimeter and tide gauge based sea level trends, and derived vertical land motion (VLM). In rightmost column compared with GPS based VLM. a) the CHON solution, b) BANH, and c) taken from Woppelmann et al. [5] Page 4 of 6

c ESA/ESTEC, 2012 20 years of Progress in Radar Altimetry, 24-26 September 2012, Venice, Italy, vol. ESA SP-710.

Figure 4. Plots of inter- and post-seismic land motions as derived from the difference between altimetry and tidal data having RMSE values ranging from low (Ko Lak and Geting), moderate (Sattahip) and large (Ko Mattaphon) Figure 4 gives the vertical land motion from the difference between altimetry and tide gauge data: blue crosses the monthly averages and the red ones 12-month moving average filtered. A linear trend has been fitted. To avoid mixing inter- and post-seismic signals, data of 6 months before and after the earthquake are not used in the fitting. Except for Ko Mattaphon, the pre-earthquake uplifts are detected at all stations and it can be seen that variations of smoothed monthly differences correlate with RMSE. The post-earthquake downward motions are also detected but the 2005–2008 time series are considered rather short to reliably define post-seismic characteristics and so the calculated trends should be regarded indicatively. The geophysical modelling shows that these large downward rates will relax over time. The contradictory behaviour of Ko Mattaphon with the other stations can be attributed to the large separation between the TOPEX-class altimetry ground track and the tide gauge (117 km), so the land motion signal has been corrupted by oceanic signals registered differently by the tide gauge and the altimeter. In Figure 5 we take a closer look at the modelling of the post and co-seismic vertical land motions. In the top panel the subsidence (post-seismicity) is given at three CGPS stations: RYNG (located near CHON), CPNT (near BANH) and ARAU (Andaman Sea, near Malaysia), both for the by GPS observed subsidence (drawn line) and the model predictions (dashed line). The subsidence rate depends mainly on the rheology of the as-

thenosphere, and over the last 7 year we see averages of about 10mm/yr. In the bottom panel we plotted the predictions for future vertical subsidence for seismic cycles of respectively 170 and 500 years by modelling the visco-elastic relaxation in the asthenosphere, as briefly explained in Section 2. Apparently, the post-seismic subsidence is temporary and will not exceed 50-60mm. Of course this is an important finding if you want to make reliable relative sea level change forecasts for certain regions. Finally we present our Interferometric SAR analyses for the greater Bangkok area in Figure 6. The plot shows the InSAR result, based on 23 RADARSAT-1 images (20052012) using StaMPS. The land subsidence rate varies between 15-25mm/yr. This resembles rates found for Venice last century and is due to water extraction and compaction. This clearly poses the largest threat in the whole Bangkok relative sea level change scenario.

4.

CONCLUSIONS

We have investigated the vertical land motion in Thailand and sea level change in the Gulf of Thailand with a focus on low-lying Bangkok. Using a multitude of data and methods that either are complementary or are used Page 5 of 6

c ESA/ESTEC, 2012 20 years of Progress in Radar Altimetry, 24-26 September 2012, Venice, Italy, vol. ESA SP-710.

–18.42

mm/yr

0.00

Figure 6. Subsidence in Bangkok from RADARSAT InSAR REFERENCES Figure 5. Post-seismic vertical relaxation: model (dashed line) comparison with GPS measured (top) and prognostications for the coming 200 year (bottom)

for verification and validation we conclude that the absolute sea level rise in the Gulf threatening Bangkok is ≈3.5mm/yr ± 0.7, but can mount to 5mm/yr near the estuaries of large rivers. Two phenomena worsen the current relative sea level change situation: tectonic subsidence from the 2004 Sumatra and subsequent earthquakes has added land downfall in the order of 10mm/yr over the last 7 years, and land subsidence from sediment compaction, which InSAR analyses pinpoint at rates up to 25mm/yr at many spots along coastal Bangkok. Good news is that the post-seismic relaxation will uplift the land after some tens of years and that the total seismic subsidence will not exceed 6cm. We also demonstrated the ability to combine satellite altimetry and GPS to replace problematic tide gauges, and propose GPS to be combined with InSAR to eliminate systematic errors.

ACKNOWLEDGEMENTS We thank the EC for funding this work in frame of the Thailand-EC cooperation facility project GEO2TECDISONG. We also thank the GEO2TECDI team, especially mentioning Anuphao Aobpaet (GISTDA), Luce Fleitout (ENS), and Luciana Fenoglio-Marc (TU Darmstadt).

[1] C. Satirapod, I. Trisirisatayawong, L. Fleitout, J.D. Garaud, and W.J.F. Simons. Vertical motions in Thailand after the 2004 Sumatra-Andaman earthquake from GPS observations and its geophysical modelling. Adv. Space Res., in press, 2012. doi: 10.1016/j.asr.2012.04.030. [2] R. Scharroo, E. Leuliette, J. Lillibridge, D. Byrne, M. Naeije, and G. Mitchum. RADS: Consistent multi-mission products. In L. Ouwehand, editor, Symp. “20 years of Progress in Radar Altimetry”, 24-26 September 2012, Venice, Italy, volume ESA SP-710. ESA/ESTEC, 2012. [3] W.J.F. Simons, A. Socquet, C. Vigny, B.A.C. Ambrosius, S.H. Abu, C. Promthong, C. Subarya, D.A. Sarsito, S. Matheussen, P. Morgan, and W. Spakman. A decade of GPS in Southeast Asia: resolving Sundaland motion and boundaries. J. Geophys. Res., 112 (B06420), 2007. doi: 10.1029/2005JB003868. [4] Itthi Trisirisatayawong, Marc Naeije, Wim Simons, and Luciana Fenoglio-Marc. Sea level change in the gulf of thailand from GPS-corrected tide gauge data and multi-satellite altimetry. Global and Planetary Change, 76:137–151, 2011. doi: 10.1016/j.gloplacha.2010.12.010. [5] G. Woppelmann, C. Letetrel, A. Santamaria, M.-N. Bouin, X. Collilieux, Z. Altamimi, S.D.P. Williams, and B. Martin Miguez. Rates of sea-level change over the past century in a geocentric reference frame. Geophys. Res. Lett., 36(L12607), 2009. doi: 10.1029/2009GL038720. Page 6 of 6