Secular changes in Earth’s shape and surface mass loading derived from combinations of reprocessed global GPS networks

The changing distribution of surface mass (oceans, atmospheric pressure, continental water storage, groundwater, lakes, snow and ice) causes detectable changes in the shape of the solid Earth, on time scales ranging from hours to millennia. Transient changes in the Earth’s shape can, regardless of cause, be readily separated from steady secular variation in surface mass loading, but other secular changes due to plate tectonics and glacial isostatic adjustment (GIA) cannot. We estimate secular station velocities from almost 11 years of high quality combined GPS position solutions (GPS weeks 1,000–1,570) submitted as part of the first international global navigation satellite system service reprocessing campaign. Individual station velocities are estimated as a linear fit, paying careful attention to outliers and offsets. We remove a suite of a priori GIA models, each with an associated set of plate tectonic Euler vectors estimated by us; the latter are shown to be insensitive to the a priori GIA model. From the coordinate time series residuals after removing the GIA models and corresponding plate tectonic velocities, we use mass-conserving continental basis functions to estimate surface mass loading including the secular term. The different GIA models lead to significant differences in the estimates of loading in selected regions. Although our loading estimates are broadly comparable with independent estimates from other satellite missions, their range highlights the need for better, more robust GIA models that incorporate 3D Earth structure and accurately represent 3D surface displacements.


Chapter 1. Introduction
The continual variation in the global distribution of surface masses (atmosphere , Tregoning and van Dam (2005b)), hydrology (Crowley et al. (2008), Tamisiea et al. (2001), van Dam et al. (2001), Dong et al. (1997)) and oceans (Chambers et al. (2004), van Dam et al. (1997))) loads and deforms the solid Earth causing changes in its shape on timescales from hours to millennia that are detectable using Global Navigation Satellite System (GNSS) data. These redistributions also introduce changes to the Earth's gravity field which cannot be directly detected using GNSS; however missions, such as the Gravity Recovery And Climate Experiment (GRACE) (Tapley et al. (2004)) can be used to measure these changes ). If the redistributions can be accurately detected and quantified then additional inferences can be made about large scale mass redistribution and the Earth's elastic response to them. These inferences will affect how the terrestrial reference frame is produced, with implications for all space geodetic missions (altimetry, Satellite Laser Ranging (SLR) and GNSS), Glacial Isostatic Adjustment (GIA) and sea level determination, amongst other applications.
Although there have been previous studies into surface mass loading , van , Wu et al. (2003)), plate tectonics , Kogan and Steblov (2008), Sella et al. (2002), Larson et al. (1997)) and GIA (Ivins and James (2005), Paulson et al. (2005), Peltier (2004), ) there is still debate as to the accuracy of some of the reference models used, for example continental water storage (CWS) (Fiedler and Doll (2007)) and GIA ), which are very difficult to model accurately.
One of the drivers behind these investigations is to gain a better understanding of the effects of climate change on the redistribution of surface masses (IPCC report 2007 (Solomon et al. (2007))).
This study utilises data from the US Global Positioning System (GPS). GPS satellites are tracked by ground stations whose position on the Earth can then be calculated on a regular basis. Regular calculation of station positions allows for the interpretation of station movements to be made. Since 1994 the International GNSS Service (IGS) has been collating data submitted by Analysis Centres (ACs) and organising processing strategies. Each week to the present, an official global network has been produced; however, clear systematic 2 processing errors are visibly present in this operational data set Ray et al. (2008).
The main issues to be addressed by this work are twofold: the first subject is the comparison of the original operational GPS solutions with the newly released reprocessed data set, detailed in section (3.6). The second focus is the level of precision and accuracy to which secular changes in the shape of the solid Earth, caused by surface mass loading, can be detected using a modern geodetic space technique such as GPS, section (6.4). Highlighted in the first issue are the implications for the Terrestrial Reference Frame (TRF), as distinct differences between the origin and orientation of the operational and reprocessed data sets have already been observed , Steigenberger et al. (2006b)). By comparing the two data sets it should be possible to highlight areas where significant geophysical signals or model deficiencies have been resolved or may still exist. To achieve this the second part of this project uses the reprocessed GPS data to develop an integrated loading model combining the secular effects of plate tectonics and GIA.
Potential error sources can be found at the very fundamental level of the data that I am proposing to use; GPS data relies heavily upon both its chosen frame of reference for the determination of its orbits and the processing model chosen.
Any error or change in the reference frame definition will filter down, contaminating the coordinates of the receiving stations, especially their heights.
Periodically a new official release of the International Terrestrial Reference Frame (ITRF) is created , Altamimi et al. (2007), Altamimi et al. (2002)); this combines information from GPS, Satellite Laser Ranging (SLR), Doppler Orbitography and Radiopositioning (DORIS) and Very Long Baseline Interferometry (VLBI). Even though its creation is based on the same principles each new release will have slight differences, in part due to changes in the observation analysis models used, section (3.3).
Over the last two decades several different TRFs and processing models were used as part of the operational GPS processing which led to an inhomogeneous coordinate time series full of systematic offsets and discontinuities, which hindered studies into long term changes in the Earth system. This was identified and in 2005 a complete reprocessing of all the GPS data was proposed by the IGS (Steigenberger et al. (2006a)) with the aim of producing a new homogeneous data set aligned to the most recent release of the ITRF, or its 3 GPS only realisation, the IGS05 (Ferland (2006)). This has since been superseded by the ITRF2008 and IGS08. The reprocessing was carried out by several ACs around the globe and part of the aim of this project is to create a new TRF from this reprocessed data series (Steigenberger et al. (2006b)). The reprocessed frame can be compared to previous frames; thus highlighting problems and model changes which have occurred with previous frames. Any issues that do arise will have wide ranging implications for many users, especially those that require the high level of precision of the ITRF.
The IGS reprocessing campaign aims to remove some of the doubts in the data sets and it is the intention of this work to develop a framework from which secular change in the Earth's shape caused by surface mass loading can be observed. Once this framework has been rigorously tested the result will be used to produce GPS derived coefficients of the Earth's deformation and converted into an equivalent mass change.
"Secular period" for geophysical studies usually refers to changes which are multi-decadal and greater, as described in Figure 1.

Aims of the Study
In this study I aim to validate and use the homogeneous IGS Reprocessing data set to develop a better understanding of surface mass loading unaffected by GPS processing model changes, e.g. the adoption of a new TRF or the switch from relative to absolute satellite antenna phase centre offsets/variations (Schmid et al. (2005)). By removing modelled GIA and estimated plate tectonic rotations, which are the two identified primary sources of secular motion, from the estimated station velocities, I will attempt to quantify the remaining secular motion caused by present day surface mass loading. In terms of this study, "secular" will be limited by the 10 years of GPS data used but long term trends will be inferred from this data.
As part of this study I have taken over the running of the Global Network Associate Analysis Centre (GNAAC) combining the operational AC solutions into a weekly combined solution. Newcastle has been part of the IGS operational combination as a GNAAC since 1995 (Davies and Blewitt (1995)).
With the advent of the reprocessing campaign Newcastle agreed to be a GNAAC again. Newcastle uses TANYA which is a bespoke piece of software developed at Newcastle by Lavallee (2000), Davies (1997) and various other authors. TANYA is a combination software which uses Least Squares (LS) to produce a unified single solution of all the individual AC weekly submissions.
For a full description of TANYA see Lavallee (2000) and Davies (1997).
Using the combined solution produced at Newcastle it will be possible to estimate a velocity for each reference frame station. These velocities result from a combination of factors including the redistribution of surface masses, plate tectonics and GIA. The amount of data available as part of the IGS reprocessing campaign begins in 1994, submitted by a subset of ACs, (4 as of NC1 combination), but I have used March 1999 to February 2010 (GPS weeks 1000 -1570) as this contains all data from all 11 participating ACs and covers over a 10 year time period. This also bridges between the reprocessing and operational data sets, at GPS week 1459.
I aim to separate out the signals in the data that are the indicators of surface loading by producing a combined GIA and plate tectonic velocity estimate. The confounding factors which have been extensively researched are GIA (Ivins and James (2005), Paulson et al. (2005), Peltier (2004), Johansson et al. (2002)) and tectonic motion , Kogan and Steblov (2008), Altamimi et 5 al. (2007), Sella et al. (2002), Larson et al. (1997)) which assigns vertical station motions to GIA and surface mass loading and horizontal motions to plate tectonics, GIA, surface mass loading and many other negligible effects which are not considered in this work.
There are several processes that occur on the surface of the Earth which can influence the weekly estimate of any chosen site: 1. The site is located upon the plate interior and moves with the tectonic plates.
2. The site is contaminated by plate boundary movements.
3. The site is located upon the plate interior, but is subject to intra-plate deformation.
4. Site contains vertical and horizontal motion due to GIA.
5. The site is insufficiently attached to the rigid plate.
6. There is a change in surface mass loading.
This is a brief summary of the main causes of a station's apparent velocity.
Alternatively the velocity of a station may be influenced by measurement errors, for example antenna/receiver changes. I aim to mitigate the effects of (1, 2, 3 and 4) as described below. I will disregard (5) as all IGS sites should have good quality monuments IGSCB (2012). Factors (1-3) can be assessed using geological observations of tectonic faults, seismology, topography, seafloor mapping and by excluding sites within the known plate boundary deformation zones. I will attempt to account for (4) by employing a GIA model. There are several GIA models in circulation (Ivins and James (2005), Schotman and Vermeersen (2005), Peltier (2004)) all of which attempt to reconstruct the geometry of the last ice age and the resulting adjustment of the solid Earth during and after glaciation. These GIA models contain potentially large scale regionally correlated errors and are likely to poorly model the magnitude of horizontal motion ). By employing different families of GIA models and comparing their output I hope to assess the effects of these errors.
Once the velocity due to GIA has been removed (4) the residual velocities should be due to plate tectonics (1-3) and any other secular loading (6). Recent publications show that plate motions can now reliably be estimated from GPS data (Kogan and Steblov (2008)) by estimating an absolute Euler pole and the magnitude of the rigid body rotation about that pole. Developing a fully integrated model comprising the modelled GIA and estimated plate velocities 6 will leave residuals caused by present day surface mass loading trends and model error.

Reference Systems and Frames
None of the above would be possible without a high quality and stable TRF.
Nearly all geodetic applications rely upon a chosen reference frame. Any reference frame is based upon a set of ideals defined by the reference frame's associated reference system.

Reference Systems
A reference system is a spatial model which fits the motion of the Earth in a theoretical space. The reference system definition contains the origin and scale of the coordinate axis and the orientation of the frame with respect to the Earth.
In the wider global geodetic community the system defined by the International Earth Rotation Service (IERS) is recognised as the reference system of choice, known as the International Terrestrial Reference System (ITRS), although there are several other Terrestrial Reference Systems (TRS) such as the ETRS and GTRS (Boucher (2001)). Focussing on the ITRS, the origin is defined as being the centre of mass (CM) of the Earth system, including the Earth and all the mass that is on the exterior of the Earth surface (atmosphere, oceans, hydrology and ice). The scale of the ITRS is defined by the SI unit of the metre.
Finally the orientation of the TRS is equatorial with the Z-axis passing through the Earth's surface at the IERS Reference Pole (IRP) which is a position equivalent to an average historic location of the North Pole. The X-axis is set to 0° Longitude (Greenwich Meridian) with the Y-axis perpendicular to the X-axis completing the right hand system. These definitions are laid down in the IERS (2004)).

Reference Frames
With the complete definition of the TRS, in this case the ITRS, it becomes possible to realise a TRF, or the ITRF from the ideals of the ITRS ). The ITRF is reference frame which is fundamental to all space geodetic applications (GPS, VLBI, SLR and DORIS). There are three stages in the realisation of the ITRF: 1. Aligning the Cartesian axes of the ITRS to the model of the Earth through a series of rotations, translations and scaling of geodetic observations.
2. Definition of the unit of length 3. Fitting an estimation of the Earth's surface to the axes.
(1) creates a physical connection between the CM of the Earth to the origin of the TRF axes; this is not a simple procedure as the CM of the Earth is constantly changing due to the dynamic nature of surface masses. Every release of the ITRF includes the location of the origin at the chosen reference epoch and its evolution in time, e.g. ITRF2005 is defined over [1995][1996][1997][1998][1999][2000][2001][2002][2003][2004][2005]. The origin of the ITRF is defined as the CM from an average of the SLR time series which for the ITRF2005 was 13 years long ).
(2) defines the scale of the frame; the SI metre is the chosen scale (McCarthy and Petit (2004)). Although this should be a repeatable measurement each technique has slightly different systematic errors and calibrations which cause slight variations to its definition. This variability means that the absolute scale of the ITRF is the weakest part of the realisation.
(3) drapes an approximated surface of the shape of the Earth over the frame and allows coordinates to be described in geographic coordinates (latitude, longitude and height (llh)) which are much easier to interpret than Cartesian coordinate (XYZ). A best approximation of the Earth is an ellipsoid defined by the semi-major (North Pole) and semi-minor axes (equatorial).

Surface Mass Loading
Surface mass loading is caused by a body of water resting upon the surface of the Earth (or as near-surface groundwater) causing deformation of the solid Earth. It is extremely important that the temporal and spatial distributions of these surface masses be quantified to be able to understand their effects on the solid Earth. Fluids on the Earth exist in three states: solid, liquid and gas. The effect of surface mass loading and redistribution affects not only the shape of the Earth but also its rotation (axis and rate) and the gravity field.
Water that is frozen can predominantly be found at high latitudes and high altitudes in glaciers and ice caps. The long-term change in the mass balance of ice stored in glaciers is a relatively stable process, which is calculated by snow fall (input) minus melt (output) which can take several years for any major secular changes (trends) to become apparent. In comparison seasonal snow 8 fall is very difficult to measure and is a weakness of hydrology models (Wouters et al. (2011)). The secular growth and retreat of glaciers causes surface loading upon vast swaths of the Earth, the resulting adjustment after deglaciation is known as GIA, section (2.4).
Most GIA models incorporate ice change until ~3-4kyr ago (Peltier (2004)); any ice change after this period is relatively small in comparison to the LGM and will not be modelled; therefore this is a potential signal in our residuals. Table 1 summarises the observed ice changes over the last ~10 years, it is changes such as this that this study hopes to be able to detect.  (Jacob et al. (2012)) Jan2003-Dec2010 Greenland up to -230±33 (Velicogna Liquid water is stored in oceans, lakes, rivers and groundwater. The amount of water stored in land can be modelled using a combination of soil moisture and hydrological models e.g. Land Dynamics (LaD) ) but these are inaccurate as global grids are interpolated from a small number of finite point records (Tregoning et al. (2009b)). Water stored in continental lakes, rivers and wetland areas accounts for less than 0.25% of the total on the surface, but the redistribution of this stored water causes noticeable signals in station time series at a seasonal and secular period (Chao et al. (2008)). Of all water mass redistribution Continental Water Storage (CWS) is the least understood due to the difficulties in making accurate recordings. Studies by Tregoning et al. (2009b), Lavallee et al. (2007), Wahr et al. (2004), Wu et al. (2003) and van  begin to formulate ideas on how the deformation due to terrestrial water is a combination of factors.
The main drivers of mass redistribution in the oceans and seas are waves, currents and tides. Waves are driven by surface winds (Tolman (2008)) introducing energy into the surface layer of the oceans. These waves have very high frequency and have a negligible influence upon loading and will not be discussed further. The ocean currents drive a massive flow of water from one area to another; these currents are caused by changes in the temperature, salinity or the wind moving small quantities of mass. However the largest 9 causes of ocean mass redistribution are the ocean tides.  (Lyard et al. (2006)) in the CM frame, (Reprocessing (2012)) to remove the effects of short-period ocean tidal loading. The choice of model is important as their coverage varies globally and some of the shallow seas are neglected in some models.
Finally water in its gaseous state mostly lies within the lowermost layer of the atmosphere known as the troposphere and contributes to atmospheric pressure. Atmospheric pressure causes a load upon the surface of the Earth and this varies with the movement of high pressure systems and short period synoptic storms passing over an area (van Dam and Wahr (1998) van Dam et al. (1997) and van Dam et al. (1994) have all studied the effects of atmospheric pressure upon station coordinates and have concluded that the effect is dependent upon a station's latitude, with the effects being greater at higher latitudes as at low latitudes the variation in atmospheric pressure is very small (van Dam and Wahr (1998)).
The quantification of the surface displacement caused by surface mass loading is achieved via Love number theory which is discussed in section (2.2).

Glacial Isostatic Adjustment
GIA is caused by a secular change in ice mass on the surface and is Movement of mantle material will occur from this induced area of high pressure to the surrounding areas of low pressure in an attempt to restore isostatic equilibrium. This movement causes a deficit below the surface and additional subsidence of the crust.
Currently the Earth is in an interglacial period which means that the large glaciers and ice sheets have receded and most of the surface of the Earth is ice free. However, the redistribution of mantle materials is not instantaneous as the mantle is highly viscous ~10 20 -10 22 Pas (Sabadini and Vermeersen (2004)). As such, GPS observations are still recording the effects of the LGM, ~26-20kyr ) caused by the slow creep of material back into its original location (Peltier (2004)). Vertical trends of up to 10mm/yr (Sella et al. (2007)) and 8mm/yr (Bevis et al. (2009)) have been measured in GPS time series at stations in close proximity to the old ice domes. There is not only vertical motion associated with GIA, generally much smaller horizontal trends (Sella et al. (2007), Johansson et al. (2002)) can be seen in affected time series. This is caused by the movement of mantle material, lithospheric flexing and the collapse of the glacial forebulge (Fjeldskaar (1994)) which is a feature of glaciation, Figure 2. GIA not only affects surface displacements; the large accumulation of ice had an associated effect upon the oceans. The amount of water in the Earth system 11 can be thought as being a constant. The growth of glaciers causes a fall in global sea levels and a reduced load upon the ocean bottom. With glacial retreat the stored water migrated back to the oceans, causing subsidence over the oceans basins due to the increased mass of water.
The sheer size of the LGM glaciers disrupts the local gravity field (Trupin (1993)) drawing the mobile fluids in the oceans towards them causing a secondary increase in loading. GIA also affects the rotation of the Earth and the location of the rotation axis. The displaced mantle moved away from the poles towards the equator slowing down the rotation rate, and increasing the Length of Day (LOD). This effect is reversed during the interglacial period as mantle material returns towards the poles ). The position of the Earth rotation axis was also disrupted. This can be seen by the slow drift of the pole towards ~80° West which happens to be the centre of the old Laurentian ice complex (Peltier (1984)).

Plate Tectonics
Plate tectonic theory, developed during the 20 th century, took great steps in unifying the study of Earth sciences. Alfred Wegener proposed his theory of continental drift in 1912 (Wegener (1912)); this outlined methods with which the continents were able to move relative to one another. This theory took steps to explain why some continental coastlines appeared to interlock for example South America and the west coast of Africa. Modern day plate tectonics arises from the combination of Wegener's theory with the later theory of seafloor spreading (Hess (1962)) which explains how new crust is created at ridges under the oceans as they drift apart. It is known that the outer part of the Earth is comprised of lithospheric plates which are ~100km thick and it is assumed that they behave rigidly in the plate interior. Tectonic plates are able to move on the surface of the sphere at the rate of 10cm/yr (Plag et al. (2006)). Where tectonic plates meet, plate boundary processes are the primary driver in the creation and destruction of the lithosphere (Condie (1997)). It is at the boundary zones that the vast majority of earthquakes of magnitude five and above on the Richter scale occur; there is a striking pattern which delineates the plate boundaries, Figure 3. There have been several attempts to calculate the velocities and rotations of the separate tectonic plates (Bullard (1965), Hess (1962)). The standout model which has been used by hundreds of studies is the NUVEL-1 Model (Demets et al. (1990)) and its updates.  (Kreemer et al. (2003) became a reliable possibility. The success of densification has been due to the improvements in portability, the cost of the GPS receivers and the immense work of the IGS densification project (Blewitt and Lavallee (1999)), section (3.1).

Issues of GIA / Tectonic Separation
The separation of GIA and tectonically induced velocities forms a large portion of the work undertaken in this thesis. Steps have been taken to identify all the causes of horizontal and vertical trends. Current methods (Klemann et al. (2008), Sella et al. (2007), Johansson et al. (2002), Cretaux et al. (1998), Argus andHeflin (1995)) assign vertical velocities to GIA and present day surface mass loading and horizontal velocities to plate tectonics and to a lesser extent GIA and surface mass loading, meaning that there is not a clear divide between the two. By ensuring that study stations are located upon the rigid plate interior it is possible to ensure that plate tectonics does not contribute to vertical velocities (Lavallee (2000)). This is further defended by Euler's theorem (Palais et al. (2009) (2004)). This is a major limiting factor to the accuracy with which the horizontal velocity can be estimated as current stratified Earth models describe each layer as having homogeneous viscosity and lithospheric thickness. If the lateral viscosity of the Earth varies ) then current models will not be able to accurately represent the observed horizontal velocity. This unaccounted for velocity will introduce a bias into any velocity estimation. Steps have been taken to adapt a stratified model by introducing areas of varying viscosity within each layer and varying the thickness of the lithosphere ), but these models are still in their infancy , Klemann et al. (2008), Latychev et al. (2005b, Kaufmann et al. (2000)). I will attempt to mitigate these GIA model shortcomings by implementing a suite of different GIA models developed around different Earth models and ice histories.

Summary
By using the newly published reprocessed IGS data set there is opportunity to study the long term deformation of the Earth. After accounting for plate tectonics and GIA, the residual velocity will be caused by present day surface mass loading trends. In the following chapters I will explain the various interpretations of GIA and how plate tectonic velocities are estimated. Chapter 3 will explain the combination of AC network solutions and modifications made to improve accuracy and inclusion of additional parameters. The final Chapters will highlight the results of the combination and how they are interpreted to quantify present day surface mass loading.
This work extends previous knowledge by using the IGS's reprocessing campaign and presents new GPS derived loading coefficients of secular mass change. By using basis functions ) and a suite of global GIA models I aim to provide several estimates of surface mass loading which can be compared against other studies, such as Wu et al. (2010), which employs an ocean model to constrain the estimation and a compilation of GIA models, or GRACE estimates (Baur et al. (2009), Velicogna (2009), Luthcke et al. (2008, Wouters et al. (2008), Chen et al. (2006b), Velicogna and Wahr (2005)).
Previously only Lavallee et al. (2007) used the basis functions but that study was concerned with extracting low degree harmonics, I am more concerned with the highest maximum level attainable and thus the finest spatial resolution.
These new results are not immune to GPS errors but they are less susceptible than the original operational data set. For example, tropospheric errors have been considered but 2 nd order ionospheric effects are not included in the first IGS reprocessing campaign. There are still GPS orbit errors which affect all GPS station positions but compared with the operational data set massive improvement is shown. With respect to loading, the homogeneous nature of reprocessing allows the study of real secular signals as there should no contribution to the long term from processing models.

Chapter 2. Geophysical Causes of Secular Deformation
Deformation of the Earth occurs over a range of time scales, from the very high frequency sub daily (tidal) to the extremely long spanning centuries to millennia.
The focus of this study is the secular signals present in GPS coordinate time series; the short term and inter-seasonal variation can be removed during the analysis by estimating a linear trend, section (4.1), and seasonal signals will average down to zero after 2.5 years of observations ), leaving only the secular deformation.
In this chapter I will discuss secular deformation; of which the two primary drivers are plate tectonics and GIA. I aim to treat these two signals as errors (biases) in the reprocessed station time series which can be quantified and removed. The linear trend of the residual coordinates, if any, will be indicative of present-day secular deformation due to surface mass loading or other processes which I will consider to be negligible. First we must consider the construction of the solid Earth and its properties as this will dictate how it reacts to surface loading.

General Earth Structure
The physical and chemical structure of the Earth has been estimated using a variety of seismological techniques (Miller et al. (2009)). By defining the layers of the Earth by their rheological properties (deformation type and strength) instead of physical and chemical properties we can begin to explain the process of plate tectonics.
The structure of the Earth can be broken down into distinct physical layers, briefly these are the crust (70km thick), mantle (2800km), outer core (2300km) and inner core (1200km). However by the rheological properties the crust and mantle can be broken down into three distinct layers (Miller et al. (2009) (2004)). The lithosphere is divided into tectonic plates; there are approximately 15 major and several other minor/micro plates ). The lithosphere has a finite strength and will rupture when placed under stress where the tectonic plates collide (Frisch et al. (2010)). Away from these rupture zones the lithosphere can be treated as elastic (Anderson (2007)) where the lithosphere is able to flex under loading. Beneath the lithosphere is a layer known as the asthenosphere, this layer is approximately 180km thick (Low Velocity Zone (LVZ) for Earthquakes) (Condie (1997)). Kearey et al. (2009) defined the root of the asthenosphere as being the beginning of the seismic transition zone. Due to the immense pressure and temperature in the asthenosphere the mantle material is ductile and it is involved in both the process of plate tectonics and GIA. Below the asthenosphere is the mesosphere, this is the portion of the mantle between the asthenosphere and the outer core (2900km depth).
Although the mesosphere is part of the mantle the increased pressure makes it much less deformable (Poirier (2000)). The Earth structure is summarised in

Elastic Deformations of the Solid Earth
The movement of surface masses on an elastic Earth will induce horizontal and vertical deformations (Farrell (1972)). The change in shape of the Earth due to the masses' gravitational and loading pressure is described by the load Love numbers (Blewitt (2003), Love (1909)).
Love number theory is used to calculate the response of the solid Earth to the disturbing forces of surface masses and luni-solar gravitation. There are two forms of Love numbers, tidal (interactions with the Sun and Moon) and loading (surface masses). The Earth's surface moves tidally by as much as 50cm per day due to the gravitational attraction of the Sun and Moon (Lennon and Baker (1973)). However this study is more concerned with the effects of surface mass loading so I will detail the load Love numbers. The (dimensionless) load Love numbers are a measure of how much a planet's surface and interior deforms in response to surface mass loading (Love (1909)). The general formula for calculating the load Love numbers is laid out in Farrell (1972). There are three different Love numbers which account for: radial displacement ' h , the additional gravitational potential ' k and the horizontal displacement ' l due to the load. If a planetary body is fully rigid then these values should all equal zero, but as the propensity of a body to deform increases then the load Love numbers increase (Love (1909)). Farrell (1972) developed a series of load Love numbers for different scenarios of Earth rheology. Blewitt (2003) building on the summary of Bomford (1980) is the unit vector in the horizontal plane and n V is the induced gravitational potential due to the surface load (Blewitt (2003)). There is an additional loading potential which is associated with the redistributed mass on and in the Earth: Using (2.1)-(2.2), the net loading potential can be shown as: This section discusses the Love numbers for an elastic Earth, however the Earth is not purely elastic and the load Love number theory can be extended to incorporate a viscoelastic Earth.

Surface Mass Loading
Surface mass loading induces deformation to the solid Earth which is recorded by GPS stations and can be described as a weekly (or daily etc.) displacement.
It has been shown that the inversion of geodetic displacement data can be used to infer surface mass loading by fitting a set of spherical harmonic coefficients (Blewitt and Clarke (2003), Wu et al. (2003)). Using the conventions from Clarke et al. (2007) it is possible to express the total load on the surface of the Earth,   There is an imbalance of GPS observations, with areas such as North America and Western Europe being saturated with observations whereas the Pacific Ocean and other areas are very sparse in observations. This imbalance causes higher degree estimation to become unstable as the required number of coefficients rapidly increases ). The standard spherical harmonic functions do not distinguish between oceans and continents and so the spherical harmonic estimate of the load will leak from the land mass to the oceans. This led to the development of an alternative method to standard 20 spherical harmonics, ), which uses a set of basis functions to constrain the loading estimation to continental land mass; this is described in section (2.3.1).

Development of Basis Functions
It has been shown that fitting a load using the standard spherical harmonic coefficients is inefficient as it suffers from the imbalance of tracking stations distribution ), as the spherical harmonic estimation is calculated evenly over the entire globe, whereas geodetic observations are mostly confined to the continental land mass. Continental land is predominantly located in the Northern Hemisphere limiting the availability of stable land in the Southern Hemisphere to locate tracking stations (e.g. GPS, SLR or DORIS).
This bias increases the uncertainties in areas which are sparse in observations or allows data sparse areas to have unrealistic values to ensure a smooth fit elsewhere.
The inverse barometer (IB) assumption Wunsch and Stammer (1997) describes how the oceans will redistribute from one geopotential surface to another by redistributing the load across the entire basin over longer time periods (greater than a few days). As such Clarke et al. (2007) strove to develop a gravitationally consistent, mass-conserving set of basis functions as an alternative to the spherical harmonic functions. These basis functions enable the estimation of realistic surface loads within the coastline, with the only signal over the oceans being the mass conserving tidal response. Clarke et al. (2007) achieved this by taking the original spherical harmonic functions and masking them with an ocean function    C which is zero over land and unity over the oceans.
As it stands (2.8) is no longer mass conserving so in the final stages of basis function development an oceanic term is introduced.
In the updated ( and is required to fulfil the SLE. The SLE accounts for the spatially varying term (change in the shape of the geoid) and a spatially constant term accounting for the total change in ocean mass.
The modified basis functions, which account for the SLE, will constrain the surface load to the continents in this case without accounting for the evolution of the coastlines. For small changes in surface mass load this evolution is negligible, in contrast to GIA modelling where coastline changes are significant.
To avoid confusion through this study I will use the term "level" to refer to the level of basis functions, i.e. n on LHS of (2.8) and "degree/order" to refer to the spherical harmonic degree and order more generally, i.e. n on RHS of (2.8). In this study the basis function are evaluated as spherical harmonics and truncated at degree 45. For a full description of the development of the basis functions see Clarke et al. (2007) who truncated them at level 30. The number of terms evaluated at a given basis function level is simply: Where N is the number of parameters and n is the basis function level. By varying the level of basis function it is possible to vary the spatial resolution of the weekly estimation to isolate specific loads, see introduction (Chapter 6), in particular the low degree harmonics.

Glacial Isostatic Adjustment
GIA is caused by response of the solid Earth to the changing surface loading caused by the growth and retreat of major ice-sheets and glaciers. With the growth of the ice-sheets come an associated fall in relative sea-levels of ~130m (Latychev et al. (2005b)). During the last major ice age, areas of North America, Fennoscandia, Antarctica and other parts of Northern Europe and South America were covered by ice, (Latychev et al. (2005b)) as illustrated in Figure 5.  GIA models attempt to determine the magnitude of deformation caused by previous glaciation (Spada et al. (2006), Ivins and James (2005), Schotman and Vermeersen (2005), Peltier (2004)). The effect of GIA can be observed through three parameters:  Changes to the Earth's gravity field (Lambert et al. (2001), Larson and van Dam (2000))  Deformation of the Earth's surface (Sella et al. (2007), Johansson et al. (2002))  Variations to the Earth's rotation ) Any one or combination of these observations can provide additional constraints upon the internal structure of the Earth and the history of previous glaciations.
Other than the direct effect of the massive loading due to glaciation there is an associated change to the localised gravity potential which varies the level and distribution of the oceans. This variation in the oceans causes an additional surface load which must be modelled (Whitehouse (2009)). By measuring rates of change of gravity and crustal deformation in previously glaciated areas it is possible to calculate constraints on the Earth's viscoelastic structure.
The uplift/subsidence of land will cause a change to the local gravity field which is detectable using absolute gravimeters; which measure the acceleration of mass in a vacuum. Repeated measurements can be used to infer the secular changes in the local gravity field, such as those caused by the change in volume/quantity of mantle mass in areas of previous glaciation (Lambert et al. (1989)). The accuracy or precision with which absolute gravity rates can be determined is 0.5μGal/yr ≈1.5mm/yr (Larson and van Dam (2000)) and more recently 0.5mm/yr (Van Camp et al. (2005)).
It is possible to measure the rate of surface displacement using space geodetic techniques such as GPS, via tracking stations attached to the solid Earth (Lidberg et al. (2010), Lidberg et al. (2007), Sella et al. (2007)) or levelling networks e.g. in Fennoscandia (Fjeldskaar et al. (2000)). Current observed rates of GIA vertical displacement have been measured at over 10mm/yr (Sella et al. (2007), Johansson et al. (2002)) in Hudson Bay and 8mm/yr in the Antarctic (Bevis et al. (2009)).
Finally, any surface deformation, redistribution of surface mass or internal structure will affect the Earth's rotation, as this induces a deformation due to centrifugal potential. For example, if a glacier forms away from the polar axis then the rotation axis would drift away from the centre of mass (as can be seen in the present day drift of the rotation pole towards Hudson bay, ~80°W (Mitrovica et al. (2005), Paulson et al. (2005)).
As this study uses GPS observations of surface deformation (with potential for future incorporation of GPS Earth Orientation Parameters (EOP)), any conclusions made will be based on the measurable surface deformation caused by GIA.

Current GIA Estimation
When completing any geodetic study it is vital that GIA is considered as either an error or a valid signal in the data. In this study I will treat GIA as a source of error which can be modelled and successfully removed from GPS observations.
To be able to accurately model the effects of GIA via surface observables, a sufficiently dense and long data set in tectonically stable areas is required. Only recently with the densification of GPS networks in Fennoscandia, see Baseline Inferences for Fennoscandian Rebound Observations, Sea-level, and Tectonics (BIFROST) (Johansson et al. (2002)) and the North American network (Sella et al. (2007)) are there data sets that can begin to constrain the models using geodetic measurements. Antarctica still suffers from poor data coverage due to the cost and inaccessibility of GNSS campaigns which is why models of Antarctica GIA are currently poorly constrained. The POLENET project (www.polenet.org) is in the process of installing new permanent GNSS receivers in Antarctica and Greenland. It will be several years before any reliable estimates of secular velocities can be made from these receivers.
Current GIA models comprise an estimated glaciation (ice) history and a modelled rheology of the Earth built from geological records. Several models have been produced over the last two decades; these models predominantly have been constrained to maximise the accuracy of vertical motions (Sella et al. (2007)). Dietrich et al. (2004) determines that GPS observations provide valuable insight into recent crustal movements in Antarctica; the results are limited however by the quality of the reference frame at the time and the limited tracking stations. Thomas et al. (2011) and Bevis et al. (2009) extends this research by using an extensive GPS network which shows that model GIA uplift is often over-estimated and that recent mass loss estimates may be systematically biased, with many being too high.

Theoretical GIA Loading
The Earth is deformed in response to a change in surface loading, the effect that a load has on the Earth is described by the load Love numbers, section (2.2) to develop a finite size time evolution of the surface loads. The total surface load is a combination of ice and ocean loads: GIA models attempt to reconstruct the historical load due to ice and oceans in a single present day estimation. Below is a short introduction of the theory; for a comprehensive explanation please refer to Whitehouse (2009). The elastic term of the same theory can be applied to present day surface mass loading which could be the driver behind any present day secular deformation, if the load change is itself secular.

Ice Load
The theoretical variation to the thickness of ice   , It can be described as: Dt is the thickness of ice at  on the Earth at the epoch t and    0 , Dt is the thickness at the reference epoch. Using this information it is possible to estimate the ice load i L : Where   i ice density. By defining the ice load within a finite boundary then: Where     spatial distribution of the load and   ft is the load time history.

Ocean Load
As the ice sheets increase and decrease in size this is accompanied by variation to global sea levels as a decreasing/increasing amount of water is stored as ice on land. At any ocean point    P for a chosen epoch t the sea level SL can be calculated by: Where  ' g r distance between the Earth's CM and the instantaneous equipotential surface and  ' i r distance between the Earth's CM and the surface of the solid Earth. By calculating the sea level then it is possible to calculate the change in sea level from a chosen reference epoch. The resulting change in sea level will change the ocean load, this can be calculated by: , Stchange in sea level. Any change in sea level (2.15) results from a change in the volume of the ocean basin, variation to the surface topography caused by vertical displacements, alterations to the geoid or changes to the amount of water contained within the ocean basins. The spatial evolution of sea level is described by the SLE (Farrell and Clark (1976)). The SLE has become a fundamental part of studying GIA and surface mass loading in recent times 27 ). In simple terms the SLE describes the evolution of the ocean level required to keep the sea surface at an equipotential, with the influx and draining of the ocean basins during the formation and collapse of glaciers.
The movement of mass from continental water storage to the ocean basins is not always instantaneous; it is highly dependent upon the surrounding topography, for example in Antarctica it can be modelled as instantaneous but continental loads in North America have many stages of intermediate storage in lakes and river, and this effect must be accounted for in a step wise manner.
Any change to the mass of water in the oceans   i mt causes a change in eustatic sea level over the entire ocean area, o A , which can be identified by a distinct spatially varying fingerprint ).
Following on from (2.11) the ocean load will change uniformly:

Construction of Models
For a GIA model to be considered fully consistent then it must account for the self-gravitation of the deformation and loads, the variation of surface ice, the distribution of the oceans and the interaction between the load and Earth's rotation. This is achieved by creating a comprehensive ice history and Earth rheology.
The ice history is compiled by studying the locations of the historical ice fronts at selected time intervals from geological records (Whitehouse (2009)) and the ice thickness are also determined from geological records (Tushingham and Peltier (1991)). The amount of ice globally is assumed from the observed global relative sea level (RSL), as the amount of ice is proportional to RSL. Ice histories are an inversion of historic relative sea level histories and therefore rely upon the quality of these records; these are supplemented by estimations implied from geological observations and rely upon inferences made from these. The maximum extent of the glaciers can be determined well from geomorphological evidence; uncertainties occur when determining ice thickness and the advance of glaciers as most evidence is destroyed with glacial retreat.

28
Our knowledge of the Earth's rheology plays an equally important role as the ice sheet. Lateral changes in the lithosphere are predominantly influenced by the motions and compositional variations of the plates ).
Vertical changes are principally due to the variation of temperature as depth increases (Mitrovica and Forte (2004)). A rheological model attempts to describe how the mantle will respond to changes in local stress and strain. It is important to accurately parameterise the interior of the Earth as this will determine the behaviour of the rocks to loading. Due to the complexity of modelling most GIA models only vary the rheology in the radial direction. Initial values are taken from an Earth model, such as the Preliminary Reference Earth Model (PREM) (Dziewonski and Anderson (1981)), which provides a comprehensive Earth model calculated from seismological studies measuring the elastic properties of the Earth.
Mantle viscosity variations are strongly linked with depth and temperature gradients and it is this heat which enables the slow creep of materials. The models described above assume that viscosity varies only with depth, however this is incorrect, there are noticeable variations in the lateral viscosity profile of the Earth, which previously have not been modelled ). The station rates interpolated from the GIA models represent the effects of glaciation as a single linear value, which is only true for the short observation period of geodetic systems. In reality the rate of adjustment is reducing from the initial unloading; there was a period of elastic adjustment initially during the unloading, after this there is an extended period of viscoelastic adjustment which rate decreases with time. The variation of mantle viscosity dictates the relaxation time of the mantle after the unloading event. The reader should also be aware of the effect of defining the depth of each layer in the Earth model, a thinner lithosphere will decrease the far-field influence of the local ice sheet, but increase the amount of localised vertical deformation (Paulson et al. (2007), Zhong et al. (2003)). Studies have shown that the assumed 120km lithospheric thickness over the British Isles is too great and that this should actually be reduced to 90km ).
The majority of trends expected due to GIA are in the vertical direction; however there is some discussion as to the magnitude of GIA in the horizontal plane (Kollo and Vermeer (2010), Teferle et al. (2009), Klemann et al. (2008, , Paulson et al. (2005)). In places the glacial forebulge/hinge line 29 horizontal velocities can exceed the vertical (Sella et al. (2007)). There is evidence that current models underestimate this effect (Lidberg et al. (2010), Sella et al. (2007), Paulson et al. (2005)).
It is well documented that the strength of the GNSS system is in the determination of horizontal positions, with vertical estimates approximately three times less precise (Leick (1995)). For 1D GIA it is the vertical estimates which have the higher accuracy in comparison to the horizontal due to the construction of the Earth model (Whitehouse (2009)). It should be possible to harness the high precision of the GNSS measurements to evaluate the accuracy of the horizontal GIA estimates and highlight areas of potential mismodelling ).
The choice of reference frame is important when defining the velocity caused by GIA; most GIA models use the CE for the definition of model geocentre.
Alternative frames such as CM are used, these different models are comparable assuming that there is no on-going surface mass exchange, (Tamisiea and Mitrovica (2011)); however this is an incorrect assumption which introduces uncertainties in model comparisons due to incorrectly matching data and reference frame differences. These errors can be mitigated by ensuring that the origins of respective models match. Tanaka et al. (2011) states that for any surface loading the CM will move towards the load and the geocentre will move in the opposite direction; the resulting displacement will move the CM away from the load. After loading, the CM initially moves away from the loaded area before moving back towards the initial load centre. It is important to have a consistent and well defined reference frame origin. For a full description of origin choice see section (2.7.2).
Although in this study I am treating GIA as a source of noise in the data there is no absolute "true" model; in an attempt to mitigate this I will use a suite of 1D GIA models and compare the results to see which models are able to best model deformation in the vertical, horizontal and a combination of both components.

Global Models Used in Thesis
I have available two families of global GIA models, Peltier's ICE5G ice history with VM2 and VM4 Earth rheologies (Peltier (2004)) and Schotman's two models based in ICE3G ice history (Schotman et al. (2008), Schotman and Vermeersen (2005)  The ice history of a GIA model describes the extent and thickness of glacial ice from the LGM until the end of the glacial period. These are provided as a regular global grid at a specified time step. Both of Peltier's models use his ICE5G ice history which covers 122-3kyr BP.
 Ice history covers 122-3kyr BP (Peltier (1996), Peltier (1994)  Maximum ice thickness 3km around Hudson Bay (Peltier (2004))  Major deglaciation ends at 3kyr BP (Peltier (2004))  121kyr-LGM set at LGM extent with only height varying (Peltier (2004)) Schotman began with the ICE3G history (Tushingham and Peltier (1991)) and made the following modifications: Although deglaciation has an absolute point in these models this does not mean that the Earth is ice free, however the deglaciation point is an indication of a halt in significant mass loss. There will be continual small changes in the amount and location of ice on the Earth but these are relatively small and therefore not included in the ice histories mentioned above. However these small changes are the focus of this study and will be discussed later.

Earth Model
The Earth model in a GIA model describes the elasticity, viscosity and internal thickness of the Earth at a regular spatial interval. Peltier (2004) provides two Earth models, VM2 and VM4, convolved with the ICE5G history and Schotman Schotman and Vermeersen (2005) provides two alternate models ( Table 2) Table 3).   There is on-going debate to the true viscosity; Schotman provides an alternative Earth structure. The alternative Earth structure varies from The variation in the lithospheric thickness will dictate the spread and depth of deformation and the increase in viscosity will delay the characteristic time scale.
Schotman notes that his changes are not warranted by any data observations but they are provided as source of comparison.
Peltier's extensive Earth models (Peltier (2004)) should provide a much more precise description of the response of the Earth to glaciation due to the greater detail of layering. These different rheologies are displayed in Figure 4.  (Peltier (2004)).
Visual examination of the pattern of horizontal velocities highlights differences between the models that are much greater than for the vertical velocities. The horizontal velocities are significantly smaller than the vertical rates but these can still produce substantial variations to velocity estimation.
Each of the figures shown in section (2.4.4) are the different models used in this study. Horizontally both of Peltier's models ( Figure 8) show a distinct secular deformation radiating away from the Hudson Bay area which is not present in either of the Schotman models ( Figure 9). These differences will be discussed in section (5.1.4).

1D versus 3D GIA Modelling
Lateral heterogeneities in the structure of the Earth result from plate tectonics, temperature gradients in the lithosphere and compositional changes in the mantle ). The main focus of 1D models has been to best-fit the vertical deformation with little regard to the horizontal (Sella et al. (2007)), yet the horizontal response is just as sensitive to the variations in the viscosity of the upper mantle (Sella et al. (2007)). It is however possible that improvement in 3D velocities will improve understanding of both the ice history and Earth models. Over the last decade, studies have been investigating the effects of laterally heterogeneous (3D) Earth models, i.e. viscosity varies with depth and location , Klemann et al. (2008), Kendall et al. (2006), Latychev et al. (2005a), Paulson et al. (2005)). It should be noted that during processing no 3D model was available for comparison but below is a discussion of the potential influences of using a 3D model instead of the standard 1D rheology.
3D rheological models are formulated using seismic shear wave velocities which can be converted via a velocity to density scaling algorithm (Latychev et al. (2005b)). These values are converted into a temperature field which it is assumed that the viscosity field is dependent upon. Latychev et al. (2005b) calculated that for North America 50% of model elements have a viscosity ~1 order of magnitude around the Lower Mantle (LM) value of 5x10 21 Pas and 95% 37 are within ~2 orders of magnitude. For the whole Upper Mantle (UM) the viscosity value is slightly more limited (factor of ~2 in log-space), but with are large variations in the upper UM (the asthenosphere).
By varying the thickness and viscosity of the different elements and introducing lateral heterogeneities different authors have been able to produce a variety of 3D Earth models (Zhong et al. (2003), Kaufmann and Wu (2002), Kaufmann et al. (2000)). The standard 1D reference model used is generally a variation on a spherically symmetric, viscoelastic Earth using the elastic structure of PREM (Dziewonski and Anderson (1981)), with an elastic lithospheric thickness of 120km, UM viscosity of ~5x10 20 Pas and LM viscosity of ~5x10 21 Pas, convolved with an ice history, e.g. ICE3G ice history (Tushingham and Peltier (1991)), this produces rates similar to those shown Figure 8.
3D models vary the thickness of the lithosphere or the viscosity of the mantle; both have their own distinct effects upon GIA rates. Latychev et al. (2005a) introduces large scale variations to the lithosphere, (220km -70km) but with an weighted average of 120km. Radial rates can be affected by as much as ±2mm/yr with horizontal motions spreading throughout the continents, perturbing rates by ~1mm/yr. Davis et al. (2008) produced four Earth model scenarios following work from Latychev et al. (2005b) convolved with the ICE5G-VM2 history (Peltier (2004)).
 1D-VM2 elastic lithosphere of 120km and radial profile from VM2 Peltier 38 model with 76800 lateral elements (dimension ~80km) and 48 radial elements (dimension ~20-110km) and a 1D model of lithospheric thickness, 120km, UM viscosity of 10 21 Pas and LM viscosity of 2x10 21 Pas convolved with ICE3G, Paulson et al. (2005) concluded that GIA rates depend most upon the viscosity of the mantle below the load. If the observations are in made in the same place as the load (e.g. GPS) then the 3D model agrees closely with the 1D model.
Localised observations in the far-field of a deglaciation area, (e.g. sea-levels along eastern seaboard) cannot be explained by the 1D as accurately as the 3D model for the Laurentian deglaciation. For truly global observations such as GRACE or 2 J , results lean towards a best fit from a 3D model but with a heavily weighted structure beneath the local load   Klemann et al. (2008) concluded that the induced motion is ~10% of the plate motion and is greater than the precision of the plate motion determination; this is summarised in Table 5.  (Wang et al. (2008), Sella et al. (2007), Steffen et al. (2006), Gasperini and Sabadini (1990), Gasperini and Sabadini (1989)) and lithospheric thickness variation determines the radial extent of an ice sheet and its influence in the farfield (Latychev et al. (2005a), Chen and Wilson (2003)). Whitehouse et al. (2006) determined that for Fennoscandia the differences between 1D and 3D model values exceed the horizontal observational uncertainties, with lateral variations in both the near-and far-field in the lithosphere and UM having a big effect upon horizontal velocities.
The change from 1D to 3D can potentially cause a 7mm/yr difference vertically and 1-2mm/yr horizontally (Latychev et al. (2009)). However it is noted that these 3D models still require study e.g. BIFROST observations still best fits a 1D model with 2-7mm/yr uplift differences for various 3D models (Steffen et al. (2006)). Davis et al. (2008) concludes that the model rate change from 1D to 3D is significant enough to warrant further investigations into lateral variations of Earth structure. Although 3D models are likely to make some difference to my results, as there were none yet available publically I have had to use the more readily accessible 1D models. Paulson et al. (2005) recommends that for localised observations in-situ to the load 1D≈3D model accuracy, otherwise a 3D model would provide more realistic values. By using a suite of 1D models it should be entirely possible to determine a good ball park estimate of the influence and removal of GIA from GPS velocities for studying present day secular mass rates. 40

Tectonics
One of the main driving forces in the Earth system is the movement of the tectonic plates. It is assumed that the tectonic plate move with a regular velocity over millions of years (DeMets et al. (2010)). To fully understand present day secular loading it is vital to precisely estimate and completely remove any secular motion due to the motion of the tectonic plates. A good understanding of the structure of the Earth is therefore required, section (2.1). Plate tectonics is concerned with the two outermost layers, the lithosphere and the asthenosphere. The lithosphere is able to move over the asthenosphere due to the viscous properties of the lower layer. However the relative movement of the lithospheric plates is hampered by frictional forces at the plate boundaries  both. This study makes the assumption that the tectonic plates are rigid (Gordon and Stein (1992)) and the interiors do not deform readily and the only seismic deformation occurs at the plate boundaries. There is some evidence of intraplate deformation. Banerjee et al. (2008) found that there is potentially 4-7mm/yr of southern motion on the Shillong plateau on the north east of the  (2000), Cretaux et al. (1998), Larson et al. (1997), Argus and Heflin (1995)). Over time, as data spans increase, more detailed models estimating a greater number of plates become available which helps to increase the overall accuracy of the solution. The second method is geological models (DeMets et al. (2010), Demets et al. (1990), Gripp and Gordon (1990)) which use geological records such as earthquakes slips (decadalcenturies), transform rates (millennial) and spreading rates (~3.16Ma) to produce an averaged estimation of the motion of the plates.
There is general agreement between geodetic models and geological models, such as NUVEL-1A, within 95% confidence. However there are uncertainties in both the geodetic observations and the geological models, i.e. omission of some minor plates or assumption that there have been no changes to plate motions over the last 3 Ma (DeMets et al. (2010), Sella et al. (2002)) which is unlikely. Sella et al. (2002) estimate that the decade of geodetic observations is representative of plate motions for the previous 10,000yrs.
The continual movement of the tectonic plates constantly changes the location of the tectonic plates in relation to one another. Any plate model attempts to describe the relative movements between plates as well as the absolute movements of individual plates. There is no method to determine the absolute orientation of the Earth's crust in inertial space as the movement of the tectonic plates is not the only motion of the Earth's crust (Steinberger and Torsvik (2008)) and therefore the absolute location of each tectonic plate is unknown; however it is possible to define an arbitrary orientation such as the No Net Rotation (NNR) model (Argus and Gordon (1991)). NNR-NUVEL-1A is an example of a NNR reference frame. NNR defines that the summation of all the plate velocities yields zero over the Earth's surface (Argus and Gordon (1991)).
By aligning geodetic models to a NNR geological model it is possible to estimate the "absolute" Euler pole of the tectonic plates. As a by-product of this study I will produce a geodetic plate velocity model which I can compare to other models, section (5.2.1).

Euler Poles
To enable any examination of plate tectonic movements it is important to be able to describe how each plate moves across the surface of the Earth. Euler's rotation theorem states "Any displacement of a rigid body, such that the geometry of the rigid body remains unchanged, is equivalent to a single rotation about a chosen axis that runs through the fixed point" (Palais et al. (2009)). The precision of estimation is improved by both increasing the number and distribution of stations on a plate (Cretaux et al. (1998)).

Estimation of Euler Poles
The determination of the absolute Euler pole location plus the angular velocity of the plate can be used to determine the plate velocity component of each station; as the distance of the station from the pole increases then the rotational velocity of the site increases. The absolute Euler pole is ): Where   A is the estimated rotation rate and e is the unit vector along the direction of the rotation axis. For any station at position r on plate A the velocity v can be calculated using the vector product: The associated speed v can also be calculated using the station's latitude  :  The motion of any plate is dependent upon the angular velocity described in ( This theory applies to estimates on a spherical Earth, however the Earth is not a sphere and the best approximation is an ellipsoid. Lavallee (2000) determined that the difference between the two for a tectonic velocity of 100mm/yr never exceeds 0.02mm/yr. This error is well below the level of noise and can therefore be disregarded.
It is a very simple process to determine the relative velocities between two plates by holding the velocity of one plate stationary and measuring the relative motion of the second plate (Cretaux et al. (1998)). This method is useful for estimating global plate circuits and fulfilling the NNR condition of the reference frame (Kreemer et al. (2006), Argus and Gordon (1991)). To calculate the relative rotation of plate A with respect to plate B we can use: Although it is very simple to estimate the relative motions of the plates, it is not possible to directly calculate the absolute Euler poles from observations so all studies rely of the accuracy of the reference frames which employ a NNR condition (e.g. ITRF2005, NNR-NUVEL-1A) , Argus and Gordon (1991)).

Forming the Model
Where i D is a  33 rotation matrix. This can be solved for p plates (Davies and Blewitt (2000)). This model only accounts for the common horizontal motion of the plate as vertical motion is assumed to be zero. There will be vertical motion input into the system caused by GIA and surface mass loading, this will not affect the estimation but it will affect the calculated unit variances.

Station Classification
The stability of the local area around the GPS tracking station is crucially important to the robustness of the station coordinate time series used in this study. It is assumed that the interior of any tectonic plate can be modelled as a rigid body (Gordon and Stein (1992)). Any station which lies within ~100km of plate boundary tends to show signs of non-rigid motion, most likely caused by the build-up of seismic strain (Kreemer et al. (2003)). Although the actual plate boundaries shown in black in Figure 11 are very narrow, several of the deformation zones shown in red spread over large areas around the plate boundaries. Kreemer et al. (2000) determined that for a globally dense and evenly distributed station network that an individual station's classification does not matter greatly when attempting to calculate a global plate circuit. To date no study has been able to implement a network of sufficient density to be able to use deformation zone stations in global circuit studies.
Even with the increased number of reprocessed stations used in this study it still falls short of the globally dense network stated by (Kreemer et al. (2000)).
Therefore it is of paramount importance that the distinction between stable and boundary sites is made so that the sites used in the plate velocity estimations are free from plate boundary discontinuities.
Although it may be possible that a boundary site has a sufficient period of clean data to precisely estimate an accurate plate rotation it cannot be included with any confidence due to post-seismic deformation or pre-seismic strain build up.
In this study once I have selected the stations which fulfil criteria of data span (>2.5 years) ) and number of observations (>104 observations) ), I will reject any sites that are likely to be contaminated by plate boundary processes (section (4.2.2)).

Expected Present-Day Mass Changes and their Observations
The focus of the discussion above has primarily been on secular changes due to GIA and tectonics and their influence upon GPS observations. However there are other geophysical causes which could be a source of secular mass change and other geodetic systems which can observe this.

Previous Studies of Surface Mass Loading
Using very low degree spherical harmonics, Blewitt et al. (2001) found that winter loading from surface fluids induces a degree-one deformation of the Earth, compressing the poles by up to 3mm and expanding the equator by approximately 1.5mm. This deformation is reversed with the changing of the seasons. Wu et al. (2002) argue however that higher degree (n≥2) signals could alias into the low degree coefficients due to the uneven network distribution.
Again using spherical harmonics but with an a priori model of atmosphere, oceans and CWS to aid ocean mass constraint Wu et al. (2003) determined that inversion of GPS displacements can accurately determine features of largescale seasonal mass variation. This method demonstrates that GPS observations agree with the models in Eurasia, North America and Australia but disagree within South America, Africa and Antarctica, Wu et al. (2003) suggested that this is due to model deficiencies and poor coverage of GPS stations which are unable to constrain local mass distribution in the inversion process. Using forward modelling, van  compared GPS displacements corrected for atmospheric pressure with a model of continental water storage. Their findings were that although the GPS observations contain 47 systematic errors, they also showed real seasonal and secular geophysical signals due to continental water loading. Wu et al. (2010) is the stand out study which attempts to predict PDMT using GRACE, Geodetic Observations (SLR, VLBI and GPS), the ECCO Ocean Bottom Pressure (OBP) model and an a priori GIA model (ICE5G with VM2 and IJ05). This is the closest alternative method to the one described in this thesis.
Wu et al's methods follow a different line of investigation in that they attempt also to estimate the GIA contribution, but ultimately aim to end up with an estimate of secular PDMT as we do. Using the above data sets, Wu et al. evidence for additional historic ice over Greenland (additional net past ice accumulation) and significant ice melt over Alaska (effects of little ice age (Larsen et al. (2004))) which are not present in the ICE5G-VM2 model; however, no physics-based GIA modelling is done to confirm this possibility and is controversial (Olaizola et al. (2011)).
In conclusion Wu et al suggests that the kinematic method allows for the simultaneous prediction of PDMT and GIA but the inversion is sensitive to the GIA and atmospheric model used as well as the GRACE data period, especially for the low degree harmonics. These results suggest that present GIA models contain large uncertainties in previously identified areas of weakness, e.g.
Antarctica, with values overstated in Greenland and significant differences over Alaska.

48
The above studies utilise observations from GPS and with the newly released IGS reprocessing campaign GPS should be able to provide new insight into present mass change dating back to the early 1990s. However there are other systems which can provide independent observations of global mass change.

SLR
Satellite Laser Ranging uses lasers to measure the distance from the surface of the Earth to orbiting satellites fitted with special reflectors. Measurements have been calculated with millimetric accuracy and have been used to determine accurate orbits and variations to the Earth's gravity field due to mass redistribution. SLR is the primary technique for determining and monitoring the geocentre and its variations. SLR has been active for over three decades now and has been used in numerous studies of the solid Earth and surface mass loading (Chen and Wilson (2008) 2001)). The SLR system provides accurate observations as most of the satellites orbit at higher altitudes, this gives them longer life and are less sensitive to high degree spherical harmonic gravity changes ).
SLR has been used to infer changes in the low degree spherical harmonic coefficients of the Earth's gravity field for many years, particularly 22 , JJ (Earth oblateness). Long term negative trends have traditionally been attributed to the mass redistribution due to GIA (Mitrovica and Peltier (1993)). Yoder et al. (1983) showed that the 2 J could be well determined from SLR observations at -2.5 to -3x10 -11 yr -1 , this value has been confirmed by more recent studies (Chen and Wilson (2003), Cox and Chao (2002), Devoti et al. (2001)). This calculated rate was consistent until 1998 when several authors noted that the trend abruptly changed to a positive value (Earth becoming less oblate) (Chen and Wilson (2003), Cox and Chao (2002), Dickey et al. (2002)) which continued until 2002 where the trend shifted back to its original rate. This became known as the 2 J anomaly and suggested a large shift in mass redistribution able to mask the signal due to GIA. Cox and Chao (2002) suggested that this anomaly was due the El-Nino/Southern Oscillation, ENSO, affecting the ocean mass, continental water storage, the hydrological cycle and snow and ice sheets. Dickey et al. (2002) estimated that the ocean could account for ~1/3 of the 2 J anomaly and suggests that the remaining redistribution could be caused by sub-polar melt.
Chen and Wilson (2003) produced a combined model of the oceans including altimetry, OBP and the ECCO, using an observation window of 1979-2002 and accounting for GIA by removing a constant trend of -2.8x10 -11 (from Cox and Chao (2002)) and concluded that the 2 J rate anomaly was triggered by the ENSO event and continued for several years due to the long relaxation time versus climatic and sea surface temperature observations. Observations of higher degree zonal harmonics 3 J and 4 J showed some small interannual variation but no large shift is present (Cox and Chao (2002)).  used SLR amongst other methods to observe the low degree harmonics, for the degree 2 harmonics, all methods showed interannual variations with good agreement. However, Benjamin et al. (2006) suggested that this apparent anomaly may be due to mismodelling of the 18.6 year solid Earth tide in the SLR analysis. Lavallée et al. (2010) showed that using a loading model accounts for the anomaly in GPS observations but still leaves a residual in SLR observations which would suggest systematic SLR errors in 2 J . It is not only long term trends of the gravity field (GIA) which affect SLR observations. Moore et al. (2005) used 30 years of SLR from Lageos and Starlette to calculate the annual and semi-annual variations (atmosphere, cryosphere, hydrosphere and oceans) of the low-degree spherical harmonics due to surface mass redistribution. The estimated values from SLR agree well with geophysical models on a 4x4 degree grid, but this agreement is not as strong when sampled on a 6x6 degree grid.

DORIS
DORIS is a Doppler system that has been included on many satellites put into orbit (Topex/Poseidon, Jason, Envisat and Cryosat) (Jayles et al. (2006)).
Unlike GPS, the receiver is on the satellite and beacons attached to the Earth re-transmit the signal which is recorded by the satellite before being downloaded. This technique began use in 1990 and there are ~50 beacons distributed globally on all major plates, with positional accuracy of ~10cm ). Cretaux et al. (1998) used observations from DORIS to determine the velocities of 45 sites of 8 tectonic plates. It was shown that results in the stable plate interior agreed with other techniques, and the distribution of DORIS sites can improve the estimation of previously poorly 50 determined plates such as Africa. These plate rates were compared with the NUVEL-1A model and had very high correlation.

GRACE
The Gravity Recovery and Climate Experiment (GRACE) is an orbiting satellite pair which measures changes to the Earth gravitation field ).
The aim of the GRACE mission is to investigate secular and seasonal gravity changes due to GIA and present day mass changes (hydrology and cryosphere). GRACE has been shown to have a geoid height accuracy of 2-3mm at a spatial resolution of ~400km (Tapley et al. (2004)). GRACE is unable to distinguish between the GIA and present day mass change so, as with GPS, GIA must be accounted for before observations can be uses for present day mass change (Tamisiea et al. (2007)). GRACE observations are provided as monthly estimations of the Earth gravity field. This provides a perfect data set for observing seasonal and secular mass variation; however, GRACE is insensitive to the degree-1 harmonics (geocentre).

Seasonal
Tests of the ability of GRACE determination of annual loading signals have been conducted by comparing values to GPS estimates. Tapley et al. (2004) demonstrated that GRACE is able to observe the seasonal variation of the Amazon watershed and separate it from the smaller surrounding watersheds. Davis et al. (2004) showed good agreement between early observations, but van Dam et al. (2007) were unable to reproduce this good agreement, which was attributed to spurious signals in the GPS. However this lack of agreement seems to have been remedied by the GPS reprocessing campaign as both Tesmer et al. (2009) and Tregoning et al. (2009b) find better agreement than van . The strength of GRACE is shown in that it is not sensitive to site specific errors which may be problematic for GPS observations.

Secular
Although hydrological signals are predominantly seasonal, some short term inter-annual changes, for example drought, will cause multi-year trends which are detectable by GRACE (Tregoning et al. (2009a)). Crowley et al. (2006) studied the water content of the Congo, and as well as the expected seasonal signal, there was an interseasonal trend of drying observed over the 4 year period equivalent to 17mm/yr. Crowley et al. (2008) conducted the same 51 investigation over the Amazon and found the expected annual signal but no evidence of drying.
One focus of this thesis is the assessment of present day secular mass loss using GPS. There have been many studies of mass loss using GRACE observation which have widely been accepted as realistic within GIA uncertainties. Over the last two decades there has been evidence of rapid acceleration of mass loss in Greenland (Velicogna and Wahr (2005)) and Antarctica (Velicogna and Wahr (2006b)). Previously these rates were determined using remote sensing, radar altimetry and ICEsat; however, these observations were limited by their temporal coverage. The GRACE mission provided a new data set with which to study glacial mass change (Chen et al.  Table 6. For a specific analysis of individual studies see section (6.5).  (2007)) and -390Gt/yr (2003+) (Meier et al. (2007)) have been estimated from GRACE, estimations of other large areas of glaciation have also been calculated, Table 7.  (1997) Matsuo and Heki (2010) The above values will be used for comparison to the estimated GPS rates presented in this study.

Terrestrial Reference Frame
TRFs are regularly realised via observations from various space geodetic techniques. To establish a fully consistent reference frame several geodetic techniques are required as no single technique can fulfil all the required parameter definitions to the highest obtainable accuracy. For example GPS is highly dependent upon the quality of the atmospheric propagation models, the satellite orbit models and the antenna phase centre models used, whereas VLBI is not affected by orbit errors and observations are calibrated with water vapour radiometers and SLR is only affected by the dry tropospheric delay and has smaller orbit modelling errors due to the nature of the satellites. These factors mean that the origin of the ITRF has traditionally been defined from SLR measurements and scale defined from VLBI and SLR, except in the formation of the ITRF2005, as the SLR scale time series was shown to be piecewise in nature ( Figure 4 of Altamimi et al. (2007)). This choice is due to the stability of the long term estimate being better than 4mm in origin and 0.5ppb in scale (Angermann et al. (2005)). This study uses the ITRF2005 ) and its GPS-only realisation (IGS05) (Ferland (2006)). The ITRF2005 presents the coordinates 0 X of selected geodetic sites at the reference epoch 0 t , (2000.0), and a corresponding linear velocity 0 X which describes the temporal evolution of the station coordinates. It is possible to calculate the station's coordinates t X at any given epoch t using: The movement of surface masses across the surface of the Earth will displace the stations' positions on seasonal and secular time periods. These displacements will manifest themselves as variations to the time series of station coordinates; which will distort the calculated weekly reference frame.
These distortions will be visible in one or more parameters of a Helmert transformation estimated between the weekly combinations and the IGS05.
Some signals that I would expect to observe are a seasonal signal in the translation along the Z-axis; forced by the seasonal mass transfer between Northern and Southern hemispheres (Kenyeres and Bruyninx (2006), ) and an annual signal in the estimated scale parameter due to seasonal hemispheric loading combined with the non-uniform distribution of tracking stations. The solution in this study and any kinematic frame is described through Cartesian station coordinates and velocities. These observations allow for the definition of geocentric distances and angles as well as inter-station parameters. For any GPS only derived reference frame the origin definition is crucial and is described in more detail below section (2.7.1).
All frames produced in this study are mapped to the IGS05 which is the GPSonly realisation of the ITRF2005. This frame has been in general use since November 2006 (GPS week 1400).
By aligning the solutions herein to the IGS05, which indirectly aligns it to the ITRF2005 ), we apply the NNR conditions (Argus and Gordon (1991), Demets et al. (1990)). This allows for the comparison of plate models published by different authors as there are all various realisations of the same TRS.
By definition the ITRF follows the principles of the ITRS and therefore the origin in should be CM, however as the ITRF does not include models of seasonal displacements this is only true over secular periods  as the 54 seasonal displacements average out. Therefore for any seasonal studies the ITRF should be viewed as having the Centre of Surface Figure (CF) as its origin and only CM in the secular (Dong et al. (2003)). The realisation of a CM origin is indirect and relies upon precise orbit modelling and processing models.

Origin Definition of TRFs
TRF's are vital to studies of geophysical processes, which has led to many studies of the importance of a TRF's definition, especially how the origin is realised (Tregoning and van Dam (2005b)). There are several different types of origin definition each of which has its benefits and drawbacks. Blewitt (2003) and Dong et al. (1997)  which is problematic (Dong et al. (2003)). Blewitt (2003) goes to great lengths to describe the three main reference frame origins, below is a brief summary of this.
The CE reference frame is fixed to the centre of the mass of the solid Earth only. However, CE is not directly accessible from geodetic observations. Initial observations (Dong et al. (1997)) describes the CE frame as a close approximation of the CF frame (see below).
CM is the centre of the mass of the solid Earth and all the surface masses in the hydrosphere. The CM is also the centre of orbit for satellites above the Earth.
Therefore, satellite geodetic observations are naturally in the CM frame but solutions rely heavily on the accurate estimation of the satellite orbits which are 55 sensitive not only to gravitational perturbations but also non-gravitational effects such as solar radiation pressure. The internal geometry is very precise; however the external solution is less precise and is highly dependent upon the orbit modelling. A CM solution is usually transformed to a previously established CM frame by ensuring No Net Translation (NNT). Any CM solution depends upon the accuracy and temporal resolution to the frame which it is tied to, e.g.
the ITRF is more suited to secular deformation but less so for seasonal loading.
Seasonal studies tend to adopt the CF (Trupin et al. (1992)), geometrically this means that the CF is defined as if the Earth's surface was covered in an infinitely dense array of points and it is the motion of these points which are taken into account. Realistically this is achieved by using a sufficiently dense global distribution of geodetic stations averaged over the globe. The origin of CF is such that the surface integral of the vector displacement field is zero CF is a suitable alternative for techniques which cannot realise a CM frame accurately, e.g. GPS as orbit modelling limiting factor and for VLBI as no satellites are used. CF≈CE (2%), CF approximates centre of mass of solid Earth (Dong et al. (1997)). A true CF origin would be difficult to realise as the surface of the Earth is constantly evolving with the creation and destruction in the plate tectonic deformation zones. This is where the Centre of Network (CN) frame is introduced, even though a sufficiently dense network is not possible a sufficiently dense global network of sites may be used to define a network origin.
Tregoning and van Dam (2005b) demonstrated via a simple diagram the effects of an uneven network on the determination of a CM, CF and CN frame origin. Dong et al. (2003) suggested that with the improvement in geodetic observations CM and CF are no longer indistinguishable.
In the CM frame the motion of any point is caused by either surface deformation motions (plate tectonics, GIA, mass redistribution) or solid Earth motions due to mass balance changes.

Geocentre Motion
Conservation of momentum dictates that the redistribution of surface mass displaces CE relative to CM , Dong et al. (1997), Trupin et al. (1992)). Both the redistribution of surface and internal masses affects the geocentre on seasonal timescales (Dong et al. (1997)), the focus is on surface mass change as internal changes are negligible on seasonal time scales.
In the CM frame every point on the surface of the Earth undergoes 2 types of motion, -Plate tectonics, fault slip, GIA, body tide deformation (surface deformation) -Solid Earth shift due to mass balance of atmosphere, oceans, ground water and internal masses Currently the ITRF only models and removes the first of these effects.
Investigations have shown that there is unlikely to be the presence of a significant secular drift between CF and CM (Dong et al. (2003)). Over periods less than 100 years GIA, sea level change and ice cover can generate an apparent secular drift between CF and CM. If there is an apparent secular drift of the geocentre then this can affect the estimation of site position which potentially causes a bias to loading estimates.

Errors in the TRF Definition
TRF scale should be repeatable for all techniques, but in reality there are several technique specific errors. GPS, for example, is highly dependent upon the ground and satellite antenna phase centre offsets/variations (Ge et al. (2005)).
The definition of the scale plays a very important role in the definition of the origin of the TRF. For example if there is an error in the scale of the GPS then this influences the ability to determine the true distance to the orbiting satellites.
A scale rate systematically changes the heights of all stations equally. Due to the uneven distribution of oceans and continents (in general) this could alias into a translation rate error, e.g. if there was a negative shift of station heights globally, the greater number of sites in the northern hemisphere could lead to an "incorrectly" estimated negative z T trend when calculating a Helmert Transformation.
In terms of loading if there is a negative trend of station heights then this would imply an increase in the amount of mass on the land/ decrease in ocean mass; this could bias the interpretation of any mass loading estimates. Morel and Willis (2005) demonstrated that 1ppb error corresponds to 10 -9 TRF scale and ~6.4mm height for stations. The shift from relative to absolute phase 57 centre antenna models was shown to produce a ~1.86ppb change to the IGS05 reference frame scale (Ferland and Piraszewski (2009) If there is an error in the TRF that the solution is mapped to then there is nothing that can be done to remedy this. An error would appear as a residual velocity and would result in a bias to one hemisphere or all stations moving uniformly up or down.
For example, previous comparisons of the Z-axis translation between ITRF2000 and ITRF2005 documented a drift of 1.8mm/yr, ); at the time this was attributed to the poor SLR network geometry, but in truth this was due to systematic errors. With the release of the ITRF2008 (Altamimi et al. (2011)) the translation rate between the ITRF2005 and ITRF2008 is 0.0mm/yr in all three components. Conclusions from these new calculations traced back to an error in the definition of the ITRF2000 origin which will only affect solutions in the operational campaign which used the ITRF2000. This would have introduced uncertainties into the solution influencing the ability to successfully remove vertical velocities due to GIA.
In theory any satellite derived TRF should be able to realise an origin centred on the centre of mass of the whole Earth as this is the focal point of the satellite orbits. This definition is highly dependent upon the satellite orbits and any mismodelling will affect the ability of the TRF to realise a true CM origin (Urschl et al. (2005)). Altamimi and Collilieux (2009) commented that the seasonal variation displayed in the Z T component of the operational GFZ AC solution was most likely linked to an orbit mismodelling influence. Further analysis , Collilieux et al. (2007)) determined that the amplitude of the seasonal variations was too large to be true geocentre motion.
There is on-going debate as to the correct method as to the best origin for the TRF attached to the Earth and which ITRF realisation best realises this origin.

58
The first argument for the ITRF2005 is presented by Collilieux and Woeppelmann (2010) who showed that network geometry can propagate errors into the Z T component by a factor of 50%. By introducing additional stations this effect was reduced to 10%. A forward model presented in this paper using tide gauge records found that the resulting reference frame showed a difference to the ITRF2005 origin of -0.44±0.22mm/yr whereas the rate to the ITRF2000 origin was 1.36±0.22mm/yr.
The argument presented by Collilieux and Woeppelmann (2010) is countered by work of  and Argus (2007) who presented the case for the ITRF2000 origin to be closer to the true CM instead of the ITRF2005.
The expected difference between the CE and CM is predicted to be no greater than 0.5mm/yr of which potential Antarctic deglaciation could explain ~0.4mm/yr of this (Argus (2007)). Argus calculated a reference frame with the origin at CE, as this is not strongly influenced by the choice of GIA model, using data from VLBI, SLR, GPS and DORIS corrected for GIA using Peltier (2004) andPeltier (1996). Predictions of CM-CE displacements indicate that Argus' method agrees with the ITRF2000 origin (-0.5mm/yr) better than the ITRF2005 origin (1.3mm/yr) as it gives a smaller CM-CE velocity. As there is no agreement over the correct realisation of the TRF origin then it can only be kept in consideration when interpreting GPS station data (Altamimi and Dermanis (2012)).

Noise
To be able to fully understand any estimates of station velocities then an accurate estimation of associated errors is required. GPS time series have been estimated over the last 10+ years and there are several sources of error which can influence the final (daily/weekly) position. Each individual station time series will have a unique error signature and especially during the operational campaign the sources of errors can evolve or change completely, e.g.
Antenna/Receiver or processing model change (Mao et al. (1999)). This should be less of a problem with the reprocessed data set as although the error sources are still present as biases, those caused by model changes will not change with time. There are two types of noise present in the GPS time series, temporally uncorrelated White Noise (WN) and temporally correlated Coloured Noise (CN). WN has equal power at all frequencies and CN has varying power, greater at lower frequencies (Santamaria-Gomez et al. (2011)).

59
Sources of noise in GPS time series include:  monument instability (King and Williams (2009))  mismodelled satellite orbits (Griffiths and Ray (2009))  reference frame errors (Argus et al. (1999))  atmospheric loading and propagation mismodelling (Tregoning and Watson (2010))  antenna phase centre model errors (Schmid et al. (2007))  local station dependent factors (King and Watson (2010))  systematic processing error (Amiri-Simkooei et al. (2007)) The noise in GPS time series can be described by the power-law process (Williams et al. (2004)), or in the time domain which has a power spectrum of: Where f is the temporal frequency, 0 P and 0 f are normalising constants and k is the spectral index. It is typical to expect naturally occurring processes to have greater power at lower frequencies and have values of k of     31 k (Williams et al. (2004), Williams (2003a) Zhang et al. (1997)).
If the noise in the GPS time series is assumed to be uncorrelated then velocity uncertainties due to WN can be approximated using the sampling rate T , providing there are a large number of observations, N : where a is the WN amplitude and T is the time series length (Santamaria-Gomez et al. (2011)). Therefore by increasing the sampling interval or length of time series, through additional uncorrelated observations, the uncertainties due to WN are significantly reduced (potentially to unrealistic values).
For time correlated noise, specifically RW the uncertainties are: with b being the amplitude of RW (Bos et al. (2008), Williams (2003a)). Thus increasing the number of position estimates (extending the data span into the time series) has very little effect upon the uncertainties and changing the sampling interval has very little effect upon the uncertainties (Santamaria-Gomez et al. (2011), Mao et al. (1999)). Bos et al. (2008) determined that for FN the uncertainties are: c is the amplitude of FN. This means that FN uncertainties are affected by the sampling rate and number of observations but not as sensitively as WN (2.27).
Noise type and amplitude can be determined using Maximum Likelihood Estimation (MLE) (e.g. Williams et al. (2004)), Spectral Index Analysis (Mao et al. (1999), Langbein and Johnson (1997) number of parameters. This method assumes that the data is temporally uncorrelated (Bos et al. (2008)), but as described here this is not the case.
Very early studies assumed that GPS time series only contained WN but this assumption was soon dismissed. Zhang et al. (1997) determined for a regional network of 10 sites that uncertainties would be 3-6 times greater when using WN+FN against WN model only. Mao et al. (1999) came to a similar conclusion except that for a globally distributed network of 23 sites the uncertainties are likely to be an order of magnitude greater for WN+FN opposed to WN only. 61 Calais (1999) determined that for a regional network in Europe that WN+FN is the best fitting model. Amiri-Simkooei et al. (2007) investigated individual station time series which all but one had 10 years of observation and a global network of 71 sites found that the best model to fit errors in GPS time series is WN+FN.
It has been noted that random walk is also present in GPS time series which has been attributed to monument instability. Johnson and Agnew (1995) found that RW error can be up to 1 2 3/ mm yr but this can be mitigated by ensuring monuments are securely anchored to bedrock, reducing RW error to 1 2 0.4 / mm yr (Johnson et al. (2000)) and therefore likely be masked by other CN sources. Williams et al. (2004) used solutions from several different analysis centres which use much denser global networks and more modern data where no attempt had been made to reduce spatially correlated noise. Again WN+FN was clearly the dominant noise model, tested against WN+RW and WN only. Williams et al. (2004) extended the statistical testing and found that global networks were shown to generally be noisier than regional networks, southern hemisphere sites tend to be noisier than those in the northern hemisphere and comparing horizontal and vertical components found that the vertical WN+FN model is 2-3 times noisier than the horizontal components.
All authors attempt to calculate a "scaling factor", a single value that can be applied to the calculated "uncorrelated" (WN) uncertainties. Bos et al. (2008) determined that the scaling factor of Mao et al. (1999)  The network contained 275 GPS sites which had 2.5+ years data span to avoid seasonal signal velocity bias ) and underestimated uncertainties due to correlated noise of short time series , Williams et al. (2004)). Santamaria-Gomez et al. (2011) found that even in the reprocessed time series a WN+CN model was superior to WN only, and further investigation into CN 62 found that FN was superior to RW. It is noted though that this does not disprove the existence of RW (station specific) only that it is either difficult to detect due to the shortness of time series or is masked by the more dominant FN in the global solution. More specifically as the spectral index and amplitude are affected by time series length then the primary noise source cannot be due to monument noise (RW) and that noise is affected predominantly by GPS processing, quality and quantity. Williams (2003b)   Hemisphere. Santamaria-Gomez et al. (2011) confirms that this is also the case for the reprocessed data. Using this theory I will calculate a realistic scaling factor for the combined solution here at Newcastle, section (4.4.2).

Summary
In this chapter I have outlined the identifiable causes of secular deformation of the solid Earth, GIA and tectonics, and how they are measured and quantified.
However there is continual debate as to the accuracy of models used especially in GIA studies. To date there have only been limited studies into the simultaneous estimation of GIA, plate tectonics and present day secular loading, e.g. Wu et al. (2010) and Wu et al. (2003). In the following chapters I will explain how the data I use in this study begins its process from raw AC input  2004versus GPS 1994. GPS is also able to study the geocentre and has greater sensitivity to the low degree harmonics, (especially n=2).

Chapter 3. Reprocessed GPS Data and Combination
Since the inception of the IGS and the completion of the GPS satellite

The IERS and IGS
There are two bodies which oversee the running of geodetic applications; these are the IERS and the IGS which are briefly described below.

The IERS
The IERS is the governing body that serves the astronomical, geodetic and geophysical communities by producing the International Celestial Reference System, (ICRS), the International Terrestrial Reference System (ITRS) and their realisations along with EOP transformations between the ICRF and ITRF, plus all the standards, constants and models required for any geodetic processing.
The ITRF produced by the IERS is indirectly accessed in this study by aligning solutions to the IGS's GPS only realisation of the ITRF (Ferland (2006)).

The IGS
The IGS is the scientific governing body of all GNSS operations; it comprises a collection of over 200 global centres that promote the exchange of permanent GNSS tracking station data to produce a variety of products. Although only one data product is explicitly used in this study, a large number of products are utilised in the steps taken to produce that product. All operations are overseen by the IGSCB which provides the recommendations and guidelines for all GNSS operations. Using various techniques these ACs, Table 9, used the raw data collected at GNSS tracking stations to produce a weekly global network based on a subset of tracking stations. These weekly solutions are deposited at a data store, Table 12. It is these weekly solutions that are combined forming the fundamental data set used in this study. Table 10 summarises the combination centres and the combination of the reprocessed data products.

Precision and Accuracy
One cause of the poor quality estimation of the station positions and velocities in previous combinations is undetected discontinuities in station time series (Williams (2003b)). Most discontinuities can be explained by changes to station hardware or Earthquakes. It is simple to account for station hardware changes as this usually appears as a single discontinuity in the time series which can be estimated as an XYZ shift, the main challenge is detecting these discontinuities.
At present there is no standard automatic process for detecting offsets; the handling of offsets is discussed in section (4.3.2). A discontinuity caused by an Earthquake cannot be handled quite as easily; after a brittle failure of the Earth,  (2006)). There are different methods for treating discontinuities. These include excluding data before the discontinuity or after dependent upon the position of the break in the time series, estimating two separate time series or in the case of post-seismic deformation excluding a section after the event or estimating a logarithmic decay function.

Station Coordinate Time Series Quality
The quality of the station coordinate time series is dependent upon several factors, including station location (visibility of sky and satellite geometry), monument stability, high quality equipment (antenna and receiver) and minimisation of multipath environment (King and Watson (2010)). In an attempt to minimise these factors the IGS guidelines ensure that station monuments are embedded into the bedrock away from any seismically active areas or local subsidence.

Reference Frame Sites
There are hundreds of sites globally that record and submit their data to the GPS data centres (Table 12). The ACs can choose any or all of these stations to be included in their solution. However when it comes to producing a reference frame the IGS has laid down some additional guidelines for station inclusion. They use the following criteria in addition to the standard IGS criteria.
This will determine whether a site will be used as a reference frame site.
 Avoidance of closely-spaced sites  Minimum 2.5 years of data  Run by a recognised geodetic institute  Part of an established regional network  Collocation to other geodetic techniques  (2000)). Ferland andPiraszewski (2009), Zerbini et al. (2007) and Ferland et al. (2000) show that by combining solutions it is possible to detect blunders and the random processing noise of each AC is averaged out. To avoid any confusion a specific naming regime for the different operational and reprocessing campaigns was agreed upon, Table 11. Other products such as ionospheric maps and tropospheric zenith path delays are not part of this first round of reprocessing. There are plans to include these and other products into future campaigns when new more robust models are available; improving the products provided for use by the wider geodetic community. Several standards were adhered to in the reprocessing campaign to ensure that all analysis centre reprocessed products were consistent, these guidelines are laid out in the position paper (Steigenberger et al. (2006a)

TANYA Combination Software
The combination of data AC reprocessed SINEX solutions presented in this thesis is achieved using the bespoke software TANYA developed and used by Newcastle University for selected IGS operations. TANYA comprises a set of executable files written in C and is controlled by a set of front end scripts written in C Shell. These control the reading of weekly network files and writing of system-dependent files through a series of logic loops and branches. These scripts enable advanced users to create and access the full power of the TANYA modules, but there is also a set of predefined scripts for beginners; see Davies (1997) for a full description.

TANYA Processing Data Flow
There are several stages in the processing scheme from the raw input files to the combined solution. Raw solution files produced by ACs are delivered to the three IGS Global Data Centres listed in Table 12. filename but a different three character extension, Table 13 . There are many file extension possibilities but only a few are of relevance here. All processing within the program requires three core extensions; in addition to these two other file extensions are used. were required by the author to correct errors or inconsistencies in the reprocessed submissions. Once the errors are corrected the corresponding TANYA block can be constructed. Each AC will impose a certain level of constraint upon its solution file; in order to be able to produce a rigorous combination these constraints must be removed, known as loosening the block.
The relevant constraint information is contained in the a priori block. The loosening of input solutions is achieved by subtracting the a priori block from the estimate block. At this time additional, arbitrary, rotational constraints can be applied to aid the inversion process (Blewitt (1997)).
Once the blocks are loosened, a Block Scaling Factor (BSF) is applied to the files' variance components; this determines the influence that each file has upon the final combined solution. The BSF is calculated from the AC unit variance: Where  2 i is the current week's BSF,   2 1 i is the following week's estimate,  is the unit variance estimate and d is the damping factor. The damping factor limits the influence that a week's variance can have in calculating the long term BSF. Davies (1997) suggests using a value of  0.2 d which Lavallee (2000) also uses and I use here. This chosen weekly damping factor only allows a 20% influence on the calculation of the following week's BSF. A varying BSF will allow for the evolution of an AC's network solution which is not always an ideal situation as this may mask any transient errors in a particular AC solution, but this procedure is required for operational processing due to the changes in processing strategies and models. For the reprocessing solutions each AC will use the same software and models throughout the entire data set and as such the variance-covariance matrices should be homogeneous. This means that it is valid to use a calculated constant BSF for each input solution.
Once the blocks have been scaled they can be combined using the rigorous least squares combination routine in TANYA. This builds an observation and parameter model for the station coordinates that appear in at least three AC solutions and is solved by inverting the standard normal equations. Included stations are put into the "g-network" solution. Any station that is rejected is put into the "a-network" solution which can be attached to the final combined solution but has no influence over the network estimation. Once completed the blocks are written in the SINEX format and are run through data quality testing.
Rejecting a station which does not appear in a minimum of three ACs for any week prevents a rogue result from an individual AC from biasing the final velocity solution.

Current Models
The main source of heterogeneity in the operational time series was identified as the various model changes implemented throughout the history of the IGS processing (Ferland (2006) (2004)) including the IGS.

Absolute Antenna Phase Centre Variation/Offsets
The model which made the biggest improvement to the precision of the global operational solution was the implementation of the absolute antenna calibrations in GPS week 1400. For short baselines with very little height change, GPS relative phase centre corrections are sufficient as the satellites in view are at a similar elevation; however this is not sufficient for global solutions (Schmid et al. (2005)). During 1996 relative GPS antenna phase centre corrections were applied by ACs to allow for a non-spherical phase response of GPS antennas (Schmid et al. (2007)). Phase centre corrections calculated consisted of an offset between the electrical and physical reference points and a variation which was dependent upon the elevation angle of the tracked satellite (Schmid et al. (2007)). The relative corrections assumed that the phase centre variation of the reference antenna (Dorne Margolin choke ring) was zero (Schmid et al. (2007)) and other assumptions meant that there were systematic errors. In November 2006 it was decided to switch to the subsequently available absolute phase centre corrections which mitigated some of the issues inherent with using the relative corrections and would therefore not bias the solutions.

IGS05 Reference Frame
The IGS05 is the GPS only realisation of the ITRF2005, produced by the IGS ). This frame is aligned to the ITRF2005 via a 14parameter Helmert transformation (Ferland (2006)) and is based on 132 sites.
This frame is chosen because it is realised via GPS observations for use with GPS time series. The reference frame to which the reprocessed solution is aligned plays an important role in the precision of station locations and the global combined network. There have been many releases of the official reference frame to which the operational solutions are aligned to; Table 14 lists all the changes of reference frame solutions each of which introduced a systematic error in the operational time series. The first round of the reprocessing campaign uses the IGS05.snx throughout.

Atmospheric Mapping Functions and Tides
The troposphere is a turbulent layer of the atmosphere which the GPS signal must pass through to reach the ground station. Of all the atmospheric layers the troposphere contains the highest water content which is highly variable. Several models have been created to map the effect of the hydrostatic and wet atmosphere path delay on the GPS signal such as the Isobaric Mapping Function (IMF) (Niell (2001)) and the Vienna Mapping Functions (VMF) (Boehm et al. (2006b), Boehm and Schuh (2004)). Of the total atmospheric path delay 90% due to dry delay is estimated fairly accurately; the remaining 10% from the wet delay is highly variable and is difficult to estimate accurately (Elgered et al. (1991)). Previously most ACs used the Niell Mapping function (NMF) (Niell (1996)) whose coefficients are empirically derived from site location and day of year. The NMF is created from data which only uses one year of northern hemisphere radiosonde profiles; this has been shown to contain latitudinal deficiencies which are at a maximum in southern hemisphere locations (Boehm et al. (2006a) (Boehm et al. (2006a)).

Higher Order Ionospheric Corrections
Published work Petrie et al. (2011), Fritsche et al. (2005 and Kedar et al. (2003) shows that higher order (2 nd and 3 rd ) ionosphere corrections produce a small systematic change in the GPS processing results, however these corrections are not considered in this round of reprocessing.

Atmospheric Tides and Non-Tidal Loading
The S 1 and S 2 atmospheric tides have no accurate model and therefore have not been included; neither were the non-tidal loading effects of the atmosphere, ocean and hydrology (Tregoning and Watson (2011), Tregoning and Watson (2010)). As there is no consensus on these effects no corrections were implemented in this round of reprocessing Penna (2011), van Dam et al. (1997)).

Receiver Bias
There are different receiver types regarding the code tracking on the L1 frequency. C1 receivers track the C/A code and P1 receivers directly track the P-code. Some receivers track P1/P2/C1, others track C1/P2. Each instrument will have a bias with respect to C1, P1 and P2. It is common to calculate the following differential biases: P1-P2 and P1-C1 (differential code bias) (Schaer (2006)). It is the second of these biases that is accounted for in the IGS reprocessing campaign Reprocessing (2012). The values for each satellite are different at the receiver, and monthly values are calculated at the CODE IGS AC and are updated each time a new satellite is launched (Steigenberger et al. (2006a)). The differential code biases are also calculated at TUM/TUD (Steigenberger et al. (2006b)) and these could potentially be used by the IGS.
For the reprocessing campaign a single mean value was used instead of a monthly variable value, this is due to the noise level of the early weeks.

Ocean Tide Loading Model
Ocean Tide Loading (OTL) results in a deformation of the Earth due to the weight of water redistributed by the ocean tides. The ocean tides are caused by the gravitational attraction of the Sun and Moon, there are multiple periodicities due to the irregular orbits of the Sun and Moon and the ocean tides can be described as a sum of all these ocean tide periodicities (Cartwright (2000)). The summation is described by an OTL model; these usually summate the main (9-11) periods. Various OTL models were used by ACs in the operational processing:  GOT00.2 (Ray (1999))  FES99 (Lefevre et al. (2002))  NAO.99b (Matsumoto et al. (2000))  CSR4.0 (Eanes and Schuler (1999))  TPXO.6.1 (Egbert et al. (1994)) However for the reprocessing campaign it was decided that the FES2004 model (Lyard et al. (2006)), corrected for CM station motion, is globally the best model available.

IGS Tracking Network
The

79
Comparing the number of stations included by each AC in the operational campaign, Figure 13 and the reprocessing campaign, Figure  Most ACs are able to balance the number of stations on the X and Y axes, however, as mentioned previously, global networks suffer from the continent rich nature of the northern hemisphere, this corresponds to the Z-axis shown by the blue lines in Figure 15. No matter what configuration of stations used there will nearly always be a bias towards the positive Z-axis.

Analysis Centre Combination Theory
There are a number of ACs committing solutions to both the operational and reprocessing campaigns. There are currently seven ACs in the operational campaign and eleven in the reprocessing, with two GNAACs. It is the role of the GNAACs to combine these network solutions into a unified rigorous solution to act as an independent comparison for the official IGS combination. Table 9 summarises each AC's contribution to the operational reprocessing campaigns.
The IGS reprocessing campaign website contains full details of each of the contributing ACs processing strategy, http://acc.igs.org/reprocess.html. It is left to the discretion of the ACs to select methods for how they carry out the reprocessing analysis, within given guidelines and recommended models.

AC Inclusion/Exclusion
Initially all of the 11 ACs that submitted reprocessed solutions were considered; however there are only 10 actually used in this study. I have chosen not to include the TIGA solution of ULR as the residuals and estimated transformation parameters suggest that there are unstated constraints that have not successfully been removed during the combination process.

Different Techniques Used
Summarised in Table 9 are the different techniques implements by the ACs as part of the reanalysis campaign. Although there were constraints on the processing models used, ACs were free to choose what software and whether to use GPS or GPS+GLONASS. Combining all the solutions will exploit the slight variations to the processing techniques to our advantage, minimising any systematic processing noise or highlighting any particular problems with individual ACs. If all solutions were calculated using exactly the same techniques and station selection criteria then any error in station coordinates or systematic processing techniques would pass unnoticed into our solution biasing any estimation. One particular AC model/technique may be superior to those used by others but this would be difficult to determine and the variation of processing techniques highlight gross errors in certain methods. ACs are free to choose what stations they include in their solution apart from the stipulated 132 IGS05 sites, the newly release IGS08 contains 232 core sites. Each centre will supplement these sites with additional sites, which still must be a recognised geodetic quality GPS sites, to their weekly solutions.

Helmert Transformation Estimation
The parameters which are fundamental to any reference frame are the origin, orientation, scale, and their evolution with time. Uncertainties in the definition of the reference frame can propagate into the precision of the station coordinates and velocity estimations. In a perfect world free from any systematic errors, the origin and scale of the reference frame will be stable without any drift compared to the truth.
Each AC solution is defined in its own unique realisation and as such they cannot be directly compared. To be able to compare the solutions, first they must be transformed to a common reference frame of choice using a Helmert transformation. A Helmert transformation is a rigid body translation, rotation, plus scale change from one reference frame station coordinates, X to another

3) is the basic Helmert transformation for a coordinate reference frame at a chosen epoch, it is comprised of three translation parameters t , a solid body rotation
R around each axis and a single scale factor . s This is applicable for station coordinates; to calculate a transformation for station velocities then (3.3) must be modified: (3.4) accounts for the rates of the different parameters where ,' XX are station velocities in the two reference frames. The rotation parameters take the form of a 3x3 matrix.
Davies and Blewitt (2000) describe how to estimate the Helmert transformation to align the solution to an existing frame. This method can also be used to remove orientation differences when calculating station residuals. Even though there is the need to estimate a rotation between the networks no analysis of these estimates is carried out as there is no relevant information to be gained.
Rotations are loosely constrained in the combination as this parameter cannot be determined from GPS measurements. All rotations are arbitrarily defined.

Combination
The combination of AC solutions are completed at selected institutions, I will focus on two of these centres. The first combination I will examine is provided by the official IGS Combination Centre (CC). The second combination results are taken from the solution produced in-house at the Newcastle GNAAC (NC1).

Combination of AC Files at the IGS
The IGS CC, during this study, was part of NRCan; here the official IGS weekly combination was computed. As with Newcastle, the IG1 combination contains the weekly estimate of selected tracking stations; it also contains information on the Earth Rotation Parameters (ERPs). There is strong agreement between the global network coordinate solutions from NC1 and IG1; therefore I will not repeat the IG1 coordinate results here. The official IG1 solution contains ERPs which are not yet included in final the NC1 solution due to numerical instabilities in some of the input AC solutions, see section (3.6.2). I will focus on the Length of Day (LOD) and the coordinates of the instantaneous rotation pole. The pole rates provide no additional information that cannot be gathered from the pole coordinates, and the Universal Time (UT) parameter is not a truly unconstrained parameter and therefore not discussed.  Figure 16 is a full estimate of the temporal evolution of the LOD and pole coordinates using the reprocessed combination data from the GNAAC IG1. In the LOD parameter there is a regular annual repeating signal, the LOD increases during the summer months where surface masses at the poles are at a minimum and there is a shortening of the day during the winter. There is a hint of a very long, perhaps inter-decadal variation but without a more substantial data span one can only speculate about this. Inspection of the coordinates of the pole reveals a regular period to the pole location; these appear to have a consistent magnitude apart from two years where there is a significant deviation from the expected range. What is known is that the redistribution of mass will alter the moment of inertia of the Earth which alters the location of the rotation axis, this is a combination of seasonal variation and GIA (Klemann and Martinec (2009)). Figure 16 shows this variation of the rotation axis is synchronised with the changing of the seasons. Figure 17 shows the XY location of the pole and its change in position since 1997. It is possible to estimate the amplitude in each of the components ( Figure 16) using a Lomb-Scargle periodogram (Press et al. (1992)). This method is used instead of the Fast Fourier Transformation (FFT) method as the FFT requires continuous data whereas the Lomb-Scargle method allows for breaks in the input data series; this is ideal for the nature of GPS data which contains some gaps due to outliers detection (Scargle (1982)).
86  Additionally there is a high frequency spike at 27.5 days which could suggest some mismodelling of sub-daily EOP. Penna et al. (2007) showed how tidal signals can alias into spurious lower frequency harmonics, with periods from 2 weeks to semi-annual. In the X and Y poles there are the two expected peaks which correspond to the annual circular motion and the Chandler Wobble of 433 days.

Combination of AC Files at Newcastle
As discussed previously the combination at Newcastle has been achieved via the bespoke software TANYA. Previously TANYA applied fiducial constraints to a sub set of core network stations, which has been shown to introduce errors into the global network estimation (Blewitt et al. (1992)). Fiducial analysis involves selecting a subset of stations and fixing their coordinates, these stations then provide a "truth" to tie the calculated network by altering the calculated subset station position so they equal the fixed values, this can distort the network and introduce large errors. In the early stages of this study steps were taken to change the processing strategy from applying fiducial constraints to a subset of core stations, to a process of estimating and applying a Helmert transformation between the weekly solution and the chosen reference frame (IGS05). This prevents any error from being introduced from applying fiducial constraints. Section (3.5.3) describes the estimation of a 14 parameter Helmert transformation between the weekly solution and the chosen reference frame.
What is not described there is the estimation of the ERP transformation. This was implemented using the techniques described by Altamimi et al. (2005): Where i is the individual solution and , RR are the rotation and rotation rate respectively,  1.002737909350795 f is the conversion factor from UT into sidereal time and  0 1 at a daily basis. UT is shown here for completeness but is not considered in the solution as it is not a freely estimable parameter by GNSS. difference in values will be very minor numerical instabilities which are inherent in the data inversion process. As a comparison of combination quality, the initial NC1 ERP's were plotted against the established routine of the IGS CC, Figure   19- Figure 20. As discussed above, due to the instability of certain input files these parameters were not included in the final combination. However as Figure 19- Figure 20 show the NC1 combined ERPs are comparable with the IG1 combination.

Combination of Operational and Reprocessed Results
The operational global network campaign has been in existence since 1994, when the tracking network was sufficiently dense to produce a consistent global network solution. The raw data collected at these tracking stations are processed by ACs; they produce all the products stated above. This study is

Block Scaling Factors
Block scaling factors are an attempt to scale the uncertainties of each AC solution by the AC's a priori estimation. During the operational processing these values were likely to change with the adoption of new processing techniques and models. However for the reprocessing campaign the estimates of a priori variances will remain relatively constant. The variation of BSF for each AC will iterate rapidly to a stable value, e.g. CO1 (Figure 21). Even after a clear slip in the estimation, most likely due to an outlier at ~110 weeks, the BSF soon iterates back to a stable level. There is a hint of a shift in the average BSF after GPS week 1459 which is the link between the reprocessed and current operational data set. This could suggest some slight systematic differences between the two "comparable" data sets.

Number of Stations in Combination
The number of stations which are common to both the solutions and reference frame are important to a certain degree, as the number of stations increases so does the degree of freedom available to estimate a Helmert transformation.   demonstrates that the X and Y axes shows a station distribution of 50% ±5% which should limit the systematic errors due to station distribution bias , Wu et al. (2002)). However, consistently the percentage of stations on the positive Z-axis is approximately 75%. Although it is not possible to get an idealised 50% split along each axis, especially the Z-axis, without thinning out several important stations, there can be at least be a fairly consistent inter-weekly distribution; which will minimise the error introduced by this uneven distribution (Collilieux and Altamimi (2009)).

Root Mean Square error and Weighted RMS
The Lavallee (2000) and Dixon (1991).   96 This improvement will also show in the combined NC1 solution.  There is therefore the strong possibility that some periods of operational processing appear better than the reprocessing. In general the WRMS of each operational reference frame period degrades over time whereas the reprocessed solution consistently reduces the WRMS, until after GPS week 1400 without any significant spikes or offsets.

Helmert Transformation Parameters
There is a high level of inter-week variability in the operational estimated translation parameters, which may mask any shifts or signals that may be present in the data. Again there is the same pattern of offsets in the scale parameter, Figure Table 16), noting a different axis scale between Figure 27 and    Table 16. This seasonal variability and lack of trend suggests real geocentre variation on a weekly basis. This is because the IGS05/ITRF2005 doesn't represent the short term geocentre only its long term evolution.

Periodic Signal Estimation
It is not always possible to detect seasonal signals in a time series or estimated transformation parameters as they may be masked by undefined offsets or high levels of noise (Williams (2003b)). This has been a problem for the operational processing and one of the aims of the reprocessing campaign is to enable study into the seasonal behaviour of GPS free from processing errors. Strong periodic signals appear in many geophysical time series and this is apparent in the reprocessed/combined GPS time series used in this study, e.g. scale in Figure   29. The main period observed in the GPS time series occurs at the annual signal with a smaller but still significant semi-annual period. These periodic signals in the GPS time series have been documented (Kenyeres and Bruyninx (2006)). Periodic signals are not an issue for the estimation of linear trends as the effect of an annual signal averages down quickly after a minimum of 2.5 years (Blewitt and Lavallee (2002) Not all signals in the GPS time series are due to geophysical causes however.
Non-loading signals have been identified in the GPS time series which is believed to be driven by systematic error combined with the repetition of the satellite geometry as it has not been found in SLR and VLBI time series , Ray et al. (2008)). The solar year has a period of 365.25 days, but the satellite geometry with respect to the Earth-Sun system repeats itself every ~351.2 days, known as the draconitic year (Ray et al. (2008)). This can manifest itself in the GPS time series producing a regular signal at the draconitic year rather than the solar year.
There has been discussion about the origin of the draconitic period (Tregoning and Watson (2010)):  Local multipath due to the satellite-geometry repeating every sidereal day, for a 24 hour sampling period the alias period is the draconitic year (King and Watson (2010)).
 Mismodelling effect in the satellite orbits.
 Errors in the a priori IERS model for the sub-daily tidal EOP variation on GPS orbits.
It is not the aim of this thesis to investigate the source of this potential error, only to note that it is present in all IGS products.

Identification of Periodic Signals
Periodic signals appear in all station time series with a varying amplitude and phase; it is possible to identify the largest of the signals in some transformation parameters or GPS time series by a quick visual inspection. The largest signals appear as a sinusoidal wave; smaller signals will be present but may not be visible due their low amplitude compared to the larger signals. This will identify signals in the time series; even those which do not appear to have significant power to influence the time series. This technique will also demonstrate the level of noise present in the GPS time series and can be used to determine the type of noise present in the data. A Lomb-Scargle expansion comprises the fundamental frequency plus any number of harmonics. The fundamental frequency in a time series can be modelled by: Any periodic signal will have an amplitude a at the given fundamental frequency f with a given phase lag  and random errors i v , in this work the Lomb-Scargle uses an oversampling factor of 10. Additional terms can be added to (3.9) as various harmonics. In the spectral domain any significant periods in the spatial time series will appear as distinct peaks in the power series. I would expect to see peaks at the annual and semi-annual periods of the solar year. Some peaks may be present at the draconitic periods and its multiples especially two, three and six (Ray et al. (2008)).
Due to the sampling frequency of the reprocessed submissions (weekly) it may not be possible to distinguish between the two periods as their periods are very close and can appear as the same signals due to the temporal resolution, Santamaria-Gomez et al. (2011) determined the resolution of distinguishable periods as 0.1cpy and Collilieux et al. (2007) determined 0.106cpy for the studies data span, although these values are not directly transferable to this study, it shows that it may not be possible to separate the draconitic and solar periods especially at the 1 st and 2 nd harmonics for ~10 years of GPS data.

Periodic Signals in Calculated Parameters
The first set of periodograms I present (Figure 30-Figure 31), will inspect the Helmert parameters of the operational and reprocessed AC solutions respectively, followed by the operational and reprocessed combined solutions ( Figure 32-Figure 33   The annual peaks in the translation components ( ,, Figure 31 could be linked to the movement of the geocentre (Moore and Wang (2003)); this geocentre motion is caused by the seasonal transportation of surface masses especially in the Z T component ). There are other 104 significant peaks (2 nd , 3 rd , 5 th and 7 th ) present in some of the translation parameters. Ray et al. (2008) examines possible harmonics in a variety of geophysical loading models and geodetic time series and concludes that for geodetic observations it is not possible to distinguish between the 1 st and 2 nd harmonics of the solar and draconitic periods (1cpy and 1.04cpy) but the higher frequency peaks are almost certainly harmonics of the draconitic as these are not present in other geodetic techniques. From the loading models Ray et al. (2008) determines that the atmosphere has a distinct annual signal but very little power sub-seasonally. A hydrological model shows potential harmonics up until 8cpy for the solar period (1cpy) but Ray et al. (2008) concludes that the loading models show strong seasonal peaks but there is little to no evidence of sub-seasonal peaks as seen in the GPS spectra. This is further examined in section (4.5).
Considering the NC1 combined operational solutions, Figure 32, there is little evidence of distinct signals in any of the Helmert Parameters apart from Y T where there is a spike at the annual period rising above the noise.

106
Through the reduced noise values there are clear spikes in x T at the 2 nd harmonic (2cpy or 2.08cpy), z T has clear spikes at the 3 rd , 5 th and 7 th harmonic of 1.04cpy and scale has a broad annual and a suggestion of the semi-annual (2cpy or 2.08cpy).
Even with the vast improvement shown in the combination of the reprocessed data versus the operational data there will always be unmodelled uncertainty in the individual station time series and in the weekly global network. Section (2.8) discusses the type and potential impact of noise on the solution and associated uncertainties.

Summary
In this section I have shown that the quality of the reprocessed solution is much better than the operational network as predicted and verified by Collilieux et al. (2011). There is still evidence of variability in some of the transformation parameters, section (3.7); this is due to the mapping of the weekly frames to a reference frame which only models the long term CM and not the weekly variability, network geometry changes and reduced residual systematic errors.
The tests described above are only a measure of the AC and combined networks quality and stability. In Chapter 4 I will discuss the steps taken to produce a kinematic velocity model including individual station error detection.

Chapter 4. Combined Time Series
In the previous chapter I discussed the combination of AC weekly solutions into a weekly unified combined network. This is the first stage of a two stage process to estimate a kinematic reference frame containing station coordinates and velocities. The second stage combines these weekly Newcastle GNAAC solutions into a single kinematic solution, section (4.1). Kinematic solutions have been used to determine crustal deformation (Kreemer et al. (2000), Lavallee and Blewitt (2000), Blewitt and Lavallee (1999)), atmospheric loading (Tregoning and van Dam (2005b), Tregoning and van Dam (2005a), van Dam et al. (1997), van Dam et al. (1994) and hydrological loading (Crowley et al. (2008, Davis et al. (2004), Dong et al. (1997). This work builds on the previous work of Lavallee (2000); in his thesis Lavallee produces a global kinematic reference frame derived from GPS observations, in addition to this a set of plate rotation vectors were calculated.
The frame was estimated using four years of global and regional GPS combinations which were combined using a free-network approach, estimating offsets and periodic signals. Lavallée uses a fractal white noise and random walk noise model when calculating his uncertainties.

Kinematic Model
The process of estimating the weekly combined solution is described below. Both X and X are  31 linear matrices of the station's Cartesian components. Davies and Blewitt (2000).
The combined data set contains weekly coordinate estimates of selected sites, The combined solution is solved by sequential Least Squares (LS) (Cross (1992)). This method stacks (sums) the normal equation components for each weekly solution i . The sequential LS solutions take the form: where i W is the associated weight matrix. The final kinematic estimate and its variance covariance matrix can be found by: As discussed previously in section (3.5.3), each weekly solution is loosely constrained to the chosen reference frame. The initial orientation and the orientation rate of the system are loosely constrained via a priori values. Each weekly solution is aligned to the chosen reference frame, in this case the IGS05 (Ferland (2006)). To be able to estimate any absolute Euler pole in the kinematic reference frame then the NNR condition must be applied. By aligning the weekly solution to a NNR reference frame each week (IGS05) it therefore applies this condition to the final kinematic solution. A complete plate model should fulfil the condition that the summation of the estimated Euler vector weighted by the plate areas should equal zero. In practice the NNR condition is only calculated over the rigid plate area and therefore it is difficult to achieve exactly zero as the estimated plate areas do not equal the absolute surface area of the Earth due to the plate boundary zones (Argus et al. (2011)).
However, Kreemer et al. (2006) determined these boundary zones to have an insignificant effect on the NNR estimation. Once solved, the LS system output is a kinematic reference frame which has an origin, scale and orientation aligned to the IGS05. It is this method that I will use to estimate station velocities by including the relevant parameters.

Study Data
The data used in this study is taken from the weekly reprocessed SINEX solutions combined here at Newcastle, NC1, into a unified weekly solution. The weekly solutions span from March 1999 (GPS week 1000) to February 2010 (GPS week 1570), which covers nearly 11 years. This is almost three times longer than the previous study of Lavallee (2000) who only uses 4.5 years. I choose to use the combined solution opposed to any individual AC solution as it is superior in precision, this is also confirmed by Collilieux et al. (2011) and Davies and Blewitt (2000). Lavallee (2000) attaches regional network solutions to the kinematic solution to densify the solution but I employ no such strategy as there are no regional networks available for the reprocessing campaign. Even without these attached regional networks the global solution in this study is substantially more dense (greater number of stations and shorter average interstation distance within each plate) than that used by Lavallee (2000)

Processing Strategy
For each week a single global combined solution is produced, with an associated rejected network comprising stations which do not appear in at least three individual solutions. Outlier detection is carried out on both of these networks as described in Davies (1997). The residual time series standard deviation is calculated for each station coordinate component in each weekly combination. Any parameter which has a residual r times greater than the standard deviation is rejected from the offending solution, and the combination is recalculated without the offending station. This is described in Table 17.  other ACs submit an incorrect solution then the correct estimate will be rejected.

In summary this method is able to detect gross outliers in an individual AC
Week 1 AC1/AC2/AC3 Week 2 AC1/AC2/AC3 Week 3 AC1/AC2/AC3 Week 4 solution, but it is unable to detect outliers if the same error is recorded by several ACs.
As discussed in section (4.1) each weekly combined solution is achieved by solving for (4.8). Following this combination each weekly combination can be combined again following (4.3)-(4.8) to produce a single kinematic solution. The LS combination method described in section (4.1), has regularly been used for combining weekly and kinematic network solutions.

Tracking network
As  The introduction of additional station coordinate observations into the velocity estimation will improve the precision of the plate estimates. There are some weekly deviations of station coordinates from the expected linear progression due to geophysical signals, modelling and systematic errors. Therefore, it is not possible to simply estimate a linear velocity for all included stations directly from the raw combined reprocessed data set. A full description of station offsets and outliers handling can be found in section (4.3). During the kinematic velocity estimation annual and semi-annual harmonic terms are estimated, however these are not removed as removal of these at this stage would influence the time series (Lavallee (2000)).
One geophysical cause of offsets in station time series are Earthquakes. Figure   35 shows the raw kinematic solution for the IGS site Conzception (conz) which 114 shows a clear offset, beside this is an example of a "clean" time series from the IGS site Potsdam (pots). There are a large number of such sites deemed suitable for use in this study ( Figure 34), which are located in the predicted plate deformation zones. It is not possible to guarantee that these time series are free from linear or non-linear site motions caused by plate boundary processes. A large amount of work and assumptions would be required to check all boundary sites for simple offsets, (which have no associated transient deformation) caused by geophysical processes. The decision was taken to prevent any contamination of station velocities from non-linear plate boundary processes, by removing plate boundary sites from any further analysis. This is done prior to velocity estimation to save the trouble of cleaning up time series when they are never used. The removal of boundary stations leaves 172 sites (red, Figure 34); these sites form the core set which will be used to calculate this study's kinematic solution. A large number of these sites are situated on the Eurasian and North American plates but there are still sufficient sites on seven other major plates to be able to estimate plate tectonic velocities (Table 18). Now that I have estimated the station coordinate time series, section (4.1), it is possible to estimate secular (linear) station velocities. The estimated linear station velocity arises from plate tectonics, GIA, secular changes in surface mass loading, and any other undetermined secular motion. In section (5.2) I discuss the estimated linear trend of plate tectonics after accounting for an a priori GIA model.

Offsets and Outlier Detection
Even with the homogeneous reprocessed time series which has had gross outliers removed from individual AC solutions, there will still be outliers and offsets due to data mismanagement (wrong station names, incorrectly recorded equipment changes…), geophysical (earthquakes, snow accumulation…) and other unexplainable causes. To be able to fully develop a global kinematic solution it is important to filter out any station outliers and offsets which have managed to pass through the data snooping routine of the weekly combination, section (4.2.1).

Outliers
Outliers are characterised as one-off deviations away from the trend of the station coordinates. In a long term solution these deviations may have only a small effect upon the overall time series, but for an individual weekly combination they can cause a large bias to the final solution. There are two separate stages of outlier detection applied during the data process leading up to a kinematic solution. The first has already been discussed in section (4.2.1).
The second data quality testing occurs when producing the kinematic solution which deals with individual station time series (along the row, Table 17 3 . The solution iterated over 80 times for this study; only when there are no more detectable outliers at the current residual level will the software reduce the  threshold and begin rejecting smaller station outliers. Figure 36 shows detected outliers circled in red, for YELL, which are clear deviations from the expected value; these observations would be removed during the iterative processing of data along with other unobservable (humanly) but still statistically significant outliers.

Offsets
Offsets have to be handled differently to outliers; they cannot be identified in the AC combination as this is only concerned with the stability of the weekly network solution and carries out no inter-week data testing. The presence of an offset can be detected in the combination as it will materialise as a degradation of a station's RMS when aligned to a reference frame; however when they occur cannot always be so readily identified. Automation is not currently possible for offsets as they will appear as an apparent deviation of a station position from its mean which does not return to the original trend the following week. This deviation effectively offsets the mean reference value and must be taken as the new mean value for detecting outliers; as if post offset data are treated as outliers large amounts of data will incorrectly be rejected. Currently there is no established automatic process to detect offsets; it must be achieved manually through visual interrogation which is subjective and very time consuming. Currently a pilot study, led at Newcastle is investigating automatic offset detection. This has begun to produce results, and is summarised in an IUGG poster ) and at (http://www.cost-es0701.geoenvi.org/working-groups/working-group-3). The process of estimating offsets effectively introduces a new station estimate into the solution with the same velocity but with an additional coordinate estimate. This "new" station is still subjected to the same criteria of data span and number of observations; in principle any number of offsets can be introduced to a station time series, but the addition of additional offsets will weaken the overall solution.
If a segment does not reach the specified criteria then it is excluded from any estimation. Once an offset has been identified then it is beneficial to determine the cause. An equipment change would cause the offset at the epoch of change and should not affect the post offset velocity and therefore the linear velocity estimation. An offset caused by seismic activity would not only comprise an instantaneous seismic deformation but potentially also an associated preseismic stress accumulation and post-seismic deformation which would appear as a transient motion to the station position (Pollitz et al. (2011)); this provides no valid data for long term studies. Earthquakes can be identified via an earthquake catalogue which categorises all major quakes, and if possible a large chunk of data must be excluded where transient motion is present. Figure   35 is a classic example of a station with clear offsets and two very different velocities either side of the offset.
In total 70 offsets have been estimated, of which all meet the data span criteria, (Appendix B). It is possible that as I am only estimating a single velocity for each station that the length of data and number of observations required does not need to be as stringent as the criteria set by Blewitt and Lavallée, but in the interest of quality I will still apply these checks. If individual velocities between offsets were estimated, to study transient velocities, then the Blewitt and Lavallee (2002) and Lavallee et al. (2006) criteria must still be applied. Even after accounting for these offsets there is still sufficient data for all 172 stations, albeit some with multiple segments.

Site Velocity Model
After several iterations of the model, I can be confident that the coordinate time series will be free from outliers. By estimating position offsets, section (4.3.2), I can be confident in the quality of data for estimating velocities using (4.1)-(4.8).
As the number of station observations increases in time so does the redundancy of the solution; this strengthens the velocity estimation and the confidence that can be put behind it. The velocity estimation is a weighted LS linear fit through all station records which minimises the error of all estimates (4.1). Figure 37 shows the estimated horizontal and vertical velocities from the reprocessed data after estimating any linear offsets, assuming Gaussian white noise and coloured noise scaling factor, section (4.4.2). A full description of the handling of uncertainties can be found in section (2.8). These velocities are the secular motion of each station estimated from the time series of station coordinates; see Appendix A for a full breakdown of station velocities. Although not the focus of this study I will briefly discuss the presence of periodic variations on the velocity estimate in section (4.5) as there are some interesting results emerging from the reprocessed data, however the focus of the thesis will remain upon the secular trends.

Discussion of Raw Site Velocities
Examination of site velocities will distinguish between the horizontal and vertical components. Observing the solution as a global model it is clear that stations which have been assigned to a tectonic plate tend to move with a common horizontal trend. There are two linear trends which will be considered as error in the time series, these are GIA, section (2.4) and plate tectonic motion, section (2.5). In Chapter 5 I will be discussing the removal of GIA trends and the estimated Euler pole/plate tectonic velocity models. These global GIA/plate models will be evaluated and tested to see which improves the goodness of fit of the global vertical trends and the effect they have on the horizontal fit. Once the horizontal and vertical velocities for each site due to GIA and plate rotation have been estimated, it is possible to remove these from the site velocity. Any remaining secular trend will be due to present day secular loading or GIA model error Chapter 6.

Estimation of Station Noise
As discussed in section (2.8), estimation of station uncertainties in the kinematic solution assuming WN will underestimate the true influence. Santamaria-Gomez et al. (2011) determined that for ~13 years of reprocessed GPS data the mean power law noise amplitude is with the mean velocity uncertainty (0.11mm/yr), a scaling factor of (13) and the RMS of uncertainties (0.13mm/yr), a scaling factor of (10.5). However these values were determined assuming that all data sets have 572 observations and were free from offsets. By repeating the calculations after excluding sites with offsets and factoring in time series length produces the FN amplitude of ~2.62mm/yr. This leads to a range of scale factor values of 131-3.56, the mean scale factor is now 25.2 and the RMS derived scaling factor is 20.2.
All of the above scale factors fall within the expected range of . By attempting to take into account individual station length and excluding offsets the scaling factor could almost double, but this could be due to a small number of stations with very short time series compared to the rest of the sample. It is for this reason that I will use a variance scaling factor of 20 when determining the uncertainties due to FN of the NC1 kinematic solution.

Periodic Site Displacements
Real periodic variations to a station's position are caused by the transportation of mass over the surface of the Earth. These signals appear in the station's time series and the magnitude and phase depends on the sites location e.g. strong hydrological signals in tropical river basins (Crowley et al. (2008)) or annual snow fall at high latitudes (Lidberg et al. (2007)).
It is not the aim of this thesis to investigate the source of this potential error, only to note that it is present to some extent in all IGS products. Due to the sampling frequency of the reprocessed submissions (weekly) it may not be possible to distinguish between the draconitic and solar annual periods, but higher harmonics should be identifiable.
Periods which are common to all sites, irrespective of location, may factor into the reprocessed station time series in addition to geophysical seasonal periods.
The examination of individual station time series will not provide any information about systematic errors in the GPS time series as each time series will be a composition of station specific and global system effects. Dong et al. (2002) estimated that approximately half of the power in GPS time series is driven by real seasonal signals, leaving the remaining 50% unaccounted for. Ray et al. (2008) stacks multiple global station time series power spectra from operational data with the aim of eliminating localised geophysical signals, highlighting common (global) higher order draconitic spikes emerging above the background noise. Ray et al. (2008) found power at 1 cycle per year (cpy) up to the 6 th harmonic, however these periods are not strictly harmonics of 1cpy. These peaks could also explained by the draconitic year of 1.040±0.008cpy and its harmonics. Ray et al. (2008) and Collilieux et al. (2007) both find no evidence of similar high periodic signals (3 rd +) in results from VLBI, SLR or loading models, suggesting that it is a problem with the GPS system opposed to a geophysical signal.
The spectra of station time series from the NC1 combination were stacked (Clarke (2012)), to enable this, a single time series was taken and all other station time series values were interpolated to this single time series; the results are plotted in Figure 38. There is clear power at the annual solar period and potentially at the semi-annual which corresponds to the global nature of the seasonal transportation of surface mass (1.0 cpy ), after which it becomes increasingly difficult to explain the higher harmonics at 1.0cpy. Using operational data Ray et al. (2008) fits an annual and semi-annual sinusoid (1.0cpy) and presents the case that it is possible to remove these two main feature without affecting the higher harmonic peaks. Through closer examination of the strong peaks (4 th harmonic UNE and 6 th harmonic NE), Ray et al. (2008) demonstrates these can be explained by using 1.04 cpy instead of 1 cpy, also there is also some suggestion of power at the 5 th , 7 th and 8 th harmonics of 1.04cpy.  124 Ray et al. (2008) and Collilieux et al. (2007) that as a similar signal has not been detected in other geodetic techniques and that the peaks tend to coincide with a period equivalent to the repeat of the GPS satellite constellation then this would suggest a systematic error in the processing of the GPS orbits. It appears that these draconitic signals remain a problem to be addressed (hopefully) in future IGS reprocessing campaigns; however, their presence should not affect secular velocity estimates which are the main focus of this thesis.

Summary
In

Chapter 5. Estimating a Combined Velocity Model
In Chapter 3 and Chapter 4 individual AC submissions were combined into weekly combined networks, which in turn were combined into station time series with estimated velocities. These linear velocities arise from a combination of factors including GIA and plate tectonics which are regarded here as a source of error. This chapter is concerned with the quantification and removal of both these sources of error. A range of GIA models will be evaluated and utilised prior to plate velocity estimation; combining these GIA and plate velocities will form the model. Figure 37 shows the estimated vertical and horizontal station velocities.
In work by Lavallee (2000) the focus was solely on the motions of the tectonic plates, so any vertical velocities were disregarded from the estimation and the horizontal velocity was assumed to be solely due to plate tectonic rotation Plate tectonics will only move sites in the horizontal plane, following Euler's rigid body rotation theory (Palais et al. (2009) By implementing a GIA model with its associated plate velocity (      0 i r t t ) in the station displacement estimation instead of a single linear velocity, and assuming that other X is zero, (5.2) suggests it is possible to interpret the residual trends as present day elastic deformation. I will then interpret this alongside models of surface mass loading.

Construction of Combined Model
Current processing techniques were designed to remove all horizontal and vertical velocities from the station time series; this process has been modified to This combined model has to be created; the first step is to compute the horizontal and vertical velocities at each station from the chosen GIA model.
This computed velocity is removed from the raw station velocities prior to the plate velocity estimation. Then combining the GIA modelling with the respective plate estimation produces the required combined model.

Null-GIA Estimation
In this section I will present the results without taking a GIA model into consideration, this will be known as the "null-GIA" hypothesis.  Figure 40 is an estimation of horizontal plate tectonic velocities at the chosen study sites; this reconstructs the method followed by Lavallee (2000). The horizontal fit of the estimated plate tectonic velocity component accounts for nearly all of the observed horizontal velocity at most sites. I will compare this to the four GIA scenarios to quantify the influence of these GIA models on the horizontal velocity estimation and the subsequent Euler pole estimation.

Tectonic Plate Designation
Following the strategy of section (2.5.5), I have designated each station to a tectonic plate, with Table 18 showing the number of plates on the rigid interior which will be used for the plate rotation estimate:

GIA Velocity Models
The null-GIA estimation, section (5.1.2), as discussed will not account for any vertical velocity in the overall estimation. To account for the GIA influence at each site I have the four GIA models previously introduced, section (2.4.4).
Each of these models will provide a different interpretation of GIA influenced station velocities and therefore produce different Euler vector and residual secular velocity results. Figure 8 and Figure 9 show the modelled GIA vertical and horizontal velocities which I will be using to make the combined loading model.
The differences between these models will allow the evaluation of the effectiveness and impact of the different GIA models within the limitations of a small set of models. It is incorrect to assume that the model that minimises the residuals is the truth. There is no conclusive "true" estimate of present day elastic deformation so any errors in a chosen GIA model may be absorbed as unknown errors in the estimation of present day surface mass loading, see Table 20 for plate estimation RMS.

Removal of Ensemble of GIA Model Combinations
The GIA models were provided as a regular grid; using bicubic interpolation between these points it is possible to estimate each station's GIA velocity component, Figure 8 and Figure 9. By forming the vector of velocities they can be subtracted from the raw station velocities.
The GIA models are only predictions of what is believed to be the isostatic response; they are assumed to be error free so there is no covariance information. To complete the calculation in TANYA a nominal covariance matrix of zeroes is implemented, but any further analysis of covariance will not be a true test of the stations' covariance. Horizontally there is very little difference to the raw observed station velocities with the addition of a GIA model. However the GIA models account for a large proportion of the observed vertical velocities. Examining the residual values there are still some vertical velocities after removal of each of the proposed models. This could either be caused by present day secular loading or by unaccounted GIA.

Plate Estimation After GIA Removal
Working with Cartesian coordinates, the rotation axis can be described by a 3x1  What Figure 42 shows is that the inclusion of the GIA model prior to plate velocity estimation makes very little differences to the overall plate estimation, with the residuals only ~1mm/yr different to the null-GIA estimation. This is a very small percentage of the overall plate velocities. There are however some distinct orientation differences between the two families of GIA model, especially over North America, the Pacific and Australia. It is these differences that the final chapter intends on exploiting when investigating present day mass trends. These small alterations are demonstrated further in (Table 20).

Absolute Euler Poles
Several studies (2000), Cretaux et al. (1998) and Larson et al. (1997) have all attempted to estimate the absolute Euler poles by aligning them to a NNR frame such as NNR-NUVEL-1A or the more recent NNR-MORVEL by applying a Helmert transformation. This allows for the comparison of absolute Euler poles estimated by different authors as all use slightly different realisations of the ITRF/IGS. Larson et al. (1997) determined that horizontal accuracies of 1.2-5mm/yr in velocity are attainable via GPS measurements. Early studies such as Larson et al. (1997) and Lavallee (2000) attempted to calculate the velocity of eight major plates using limited GPS data with an uneven distribution. More recently there have been 133 studies which have had a much larger data set at their disposal. Sella et al. (2002) presented the REVEL model which describes the relative velocities of 19 plates, Sella et al. (2002) used operational GPS data from between 1993-2000, (reprocessed at the time to limit inconsistencies) using 200 GPS stations.
Instead of applying a particular GIA model, Sella et al. (2002) removed any sites deemed to be affected by GIA, this is achieved using the ICE4G model (Peltier (1994)) and rejecting any sites with modelled horizontal rates of ≥0.7mm/yr and vertical rates of ≥3mm/yr. Sites are rejected in North America, Eurasia and Antarctica. Even more recent is the NNR-MORVEL56 model Argus et al. (2011) which defines the angular velocities of 56 plates in a NNR reference frame. This model comprises of 25 plates from the MORVEL model (DeMets et al. (2010)) and 31 minor/micro plates from Bird (2003). The addition of Bird's 31 plates increases the coverage of the Earth's surface by 2.8%.
Comparing this frame to the results presented in this study enables the comparison of two different NNR realisations, ITRF2005 (this study) and Argus' own frame. Altamimi et al. (2007) estimates 15 major plates from a velocity field of 152 sites which have horizontal velocity errors of less than 1.5mm/yr and are located on the rigid plate interior and have 3+ years of data. Table 19 lists the common estimated absolute Euler poles plus associated error ellipses. Each of these studies is supposed to meet the NNR condition; however each of these solutions is aligned to a different reference frame. Although these should fulfil the NNR, there will be slight variations to the reference frame realisations. There are large improvements to the formal errors for the Nubian, South America and Pacific plates which are down to the increased number of stations assigned to each plate from early studies (Lavallee (2000), Larson et al. (1997)) to the more recent studies , Sella et al. (2002)). The size of the error ellipse depends on both the error estimation from each input AC but also the noise model used in the kinematic estimation,    Table 20 lists the absolute Euler poles of the nine plates listed above as in Table 19, but after accounting for the various GIA models (stars in Figure 43     suggest that the presentday GIA motion is somewhere in between the ICE4G (not shown) and ICE5G models and that the VM2 model is a poor horizontal fit to recent observations in particular over North America which would be one explanation for the values in Table 20.
Vertically the RMS clearly shows the influence that introducing a GIA model has on the estimation. First focussing on areas of known previous glaciation there is a reduction of RMS for all models in North America (15%-44%) and Eurasia (15%-25%), this greater reduction over North America is most likely due to the higher percentage of North America covered by the Laurentia ice sheet compared to the more limited Fennoscandian ice sheets coverage of Eurasia.
There is some disagreement over Antarctica (-22% to +45%) with no pattern to GIA family, this highlights the need to improve the understanding of GIA over Antarctica (Thomas et al. (2011)). Of the plates with "no" GIA the majority unsurprisingly show little to no change to the vertical RMS with the addition of a GIA models. This is true for Australia (7%), India (5%), Nubia (9%) and Pacific (15%, very small initial values), though there is a distinct divide between the vertical RMS for Nazca (+23%) and South America (-25%) for the Peltier models (-5% and -6% for Schotman). This small variation suggests some mismodelling of Peltier's GIA models either over South America (Chilean glaciation) or North America affecting the far field (King et al. (2012)).

Secular Velocity Model Combination
After estimating the GIA adjusted plate tectonic velocities the two components of the secular velocity model are now in place, one being the given GIA model and the second the associated estimated tectonic plate model. These provide velocities for each of the study sites which can be summed (5.4) to produce the required model.
The estimated horizontal GIA velocities are never more than ~15% of the corresponding plate velocities so the removal of these should only make slight variations to the residual horizontal velocities. However the residual velocities are also ~15% of the estimate plate velocities, suggesting that the plate residuals could either be due to present day surface mass loading, GIA model error or a combination of the two.

Visual Inspection of Models
The residual velocities will highlight expected areas of present day surface mass loading and/or model errors ( have an effect in Fennoscandia to a lesser degree as there is a higher ratio of land to water which is able to contain a surface load and hence a greater probability that terrestrial water storage could explain these residuals instead of model errors but due to the limited number of stations in the area is difficult to quantify.

Statistical Testing of Models
A  2 dof test will attempt to surmise the suitability of each model produced in   (Table   20). However, the alternative Schotman model also happens to have the least beneficial effect on the vertical  2 dof fit, (23%). It should be noted that even though Peltier's models significantly increase the horizontal  2 dof they will not be discarded as it is possible that the null-GIA method is absorbing some of the horizontal GIA velocity into the tectonic velocity or there is an error in the model over North America and retaining them will provide an interesting comparison of potential secular mass loading.

Summary
In this chapter I have presented the processes used to produce the combined model of estimated secular velocities due to GIA and plate tectonics. The residual velocities to each of these models are the potential result of present day secular loading. The introduction of Schotman's primary model gives the best improvement over the null-GIA estimate, however all models will still be considered. In Chapter 6 I will present an estimation of loading by attempting to fit a set of modified spherical harmonics and evaluate the suitability of each model for quantifying the calculated mass change of selected areas.

Chapter 6. Secular Loading
The previous chapter describes the construction of the surface velocity model (GIA and plate rotation). Once this modelled velocity has been removed from the raw observations, the residual velocities are assumed to result from present day secular surface mass loading. To fully understand and interpret the residual velocities, a modified spherical harmonic estimation, section (2.3.1), of the secular surface mass loading will be produced. As the degree of estimation increases, the spatial resolution becomes finer which allows for the detection of smaller-scale loading signals. As a rough guide, it is possible to estimate the spatial scale of each basis function level n .  To avoid confusion, the terms "degree" and "order" refer to the standard spherical harmonic degree and order. The term "level" refers to the spherical harmonic degree from which the alternative basis functions are derived ). The estimated basis function coefficients are reconstituted into spherical harmonic coefficients (to maximum degree and order 45, i.e. a spatial resolution of 440 km). The level of estimation can be increased until ultimately there are too many parameters and not enough independent observations to calculate a solution. It is important to ascertain the maximum level estimable from the available data sets as the distribution of stations, able to constrain the inversion, is uneven; which means the effective maximum level for each region is not the same and therefore the maximum global level must account for this.
There is no definitive method to determine the correct level to estimate a 147 loading fit. Methods include comparing the spatially re-constructed estimated load with a "true" model of the predicted surface loads, but in practice there is no way to determine the true value of loading. A second method would examine the smoothness of fitted coefficients week to week, but again this is not a rigorous quantitative measure of truth, as it is clear that surface loads cannot move instantaneously from one location to another. Thirdly, as described below, it would be possible to examine the goodness of fit of a particular estimate to the observational data. A final method, and the one used ultimately in this thesis, would be to create a synthetic data set from a realistic loading model, and test the ability of the fitted coefficients to recover the known load.

Goodness of Fit of Basis Functions
One  Once the  2 dof for each model at each level has been calculated then it will be possible to determine the optimum global level estimable. The above table is calculated using all the observations; this global statistic however implicitly assumes that the stations are evenly distributed over the surface of the Earth, which is not the case. This test shows that the increasing level of basis function globally fits the data roughly equally well for each maximum level, and within the expected range of goodness of fit, (Figure 48).
The introduction of a GIA model is having a quantifiable effect at all levels; however, it is not possible to reject any particular basis function level as an implausibly good or bad fit globally. Realistically the global dof will not represent the true level of constraint imposed by the data over the different regions. Over-fitting in one region may be 149 compensated by under-fitting elsewhere. Section (6.2) shows how this may happen.

Synthetic Data Analysis
To test the effectiveness of the basis functions, a realistic synthetic data set has been created. The data set has been created from the combined loading model of atmosphere, National Centre for Environmental Prediction (NCEP) reanalysis (Kalnay et al. (1996)), the Land Dynamics (LaD) continental water storage ) and OBP from the Estimating the Circulation and Climate of the Ocean (ECCO) general circulation model (http://www.ecco-group.org).  This total seasonal load takes the form of an annual sine and cosine grid ( Figure 49). To create a realistic loading displacement at a station then the data must be converted to a total load    T , equation (6.3), then to an equivalent height and lateral displacement using the load Love numbers, section (2.2),. Where  is the geographical location and cos sin , TT are the cosine and sine coefficients of the load. (6.3) also interpolates the UNE displacement value for each station used in this study for a single year which is repeated to produce a decadal time series. This time series is converted into a corresponding TANYA block (GPS week 1000-1570), equivalent to the loading displacements 151 produced as part of the kinematic estimate, section (4.4). The purpose of this synthetic data set is to test the sensitivity of the basis function routine to white noise and velocity biases.
This "true" data set is corrupted by adding synthesised Gaussian white noise to each weekly block, with magnitude given by each station's scaled formal position error in that weekly block. A second simulation has been compiled which additionally introduces a velocity bias; this has been achieved by producing a random set of velocity biases which have a mean of zero and standard deviations of 0.2mm/yr and 1mm/yr in the horizontal and vertical respectively. The actual range of random velocity error in each of the simulations is ±0.6mm/yr in the horizontal and ±3mm/yr in the vertical.
The data set described above and the tests below are an attempt to produce some indication of the inversion error of the basis functions. However the reader should be made aware of a few limitations of the synthetic tests.
It must be noted that the synthetic data set (Clarke et al. (2005)) does not contain any substantial mass change over the Antarctic region as the hydrological model omits it and the atmospheric component is very small. This lack of model data will hinder the ability of the inversion test to calculate a realistic error budget.
The data set is interpolated over the full network of 172 stations, in reality though the trends are calculated over a network which is constantly changing its size and distribution. This lack of temporal network variation will likely cause an underestimation of the true level of inversion artefacts. In an attempt to quantify the effect of this a reduced network was created (section 6.2.1). This reduced network size is an attempted recreation of the early GPS network to determine the effects on the final inversion. This however is only one test, in reality the very early GPS network was much smaller than the two tested here and latter weeks were much larger ( Figure 14). The reader should be made aware that the estimated error budget calculated will underestimated the true error due to network evolution.
The tests below only consider random errors, but in reality it is highly likely that there is strong spatial correlation to these errors such as those due to GIA ). These small errors with a regional correlation can have a significant influence upon the final mass estimate, by comparing the Peltier and Schotman models over North America which has a distinctly different spatial pattern the reader can begin to understand the effect that small regionally correlated errors can have on the final solution.

White Noise
For the weekly white noise simulation there should be no secular trend present at any station's loading displacements, therefore the basis function estimation should also contain no trends in any of the reconstituted spherical harmonic coefficients.
One simulation would not be statistically sufficient, several (four) syntheses were compiled, inverted and the RMS of each level's (4-6) trend was calculated,

154
Small trends of ±3mm/yr exist in the reconstituted loads especially at higher levels and in areas with fewer GPS stations. These trends can be thought of as the precision or numerical error of the basis function estimations.

Network Distribution
The conclusions stated above are for a particular network geometry using the same synthetic load. Network geometry has been shown to be a potential source of error when estimating a loading inversion ), especially the bias of stations towards the Northern Hemisphere. To test this hypothesis, a thinned network of study sites has been produced; this was achieved using Delaunay triangulation to maximise inter-site distance.
Specifically, the inter-site distances for the original network were calculated and if a site appeared within 200km of another site it was flagged. The site with the most vectors less than this threshold was then removed. This procedure was iterated until no vectors fell below the threshold (some discretion was used over southern hemisphere sites), at which point the threshold was increased by 100 km and the procedure repeated until a network of ~80 sites remained (at a final threshold of 500 km). The new network maintains a global distribution but evens out the distribution of sites between continents and especially along the Z-axis.
Following the same methods as above a set of plots ( Figure 51) was created which represents the RMS of four synthetic data sets with white noise, these solutions show more regions of uncertainties but again these are very small and disperse. Again, it would be expected that there should be no trends anywhere, but as shown by Figure 51, for a thinned out network the basis function estimation accumulates in areas which are poor in spatial coverage (e.g. Madagascar, Himalayas, no sites). However, the maximum error for all tested levels never reaches greater than ~3mm/yr globally, apart from Madagascar.

Velocity Bias
For the data set which has had a velocity bias added then each station will have a random trend added; these trends will be recovered by the basis function estimation and appear as (erroneous) trends in the reconstituted spherical harmonic coefficients. Again several (four) syntheses were compiled, inverted and the RMS of each level's (4-7) trend was calculated, Figure 52-Figure 54.
Again a contour has been chosen, this time 40mm/yr, to highlight the regions of confidence. This value is also subjective to reflect the potential divide between real and artefacts. Figure

160
The conclusions to be taken from the synthetic data set for the networks and geometry shown are:  Level 5 is the maximum global level estimable without introducing major artefacts, however areas with a good network coverage should be reliable for estimates to a higher level (7)  The areas susceptible to basis function estimation are those that are lacking tracking stations (or have had them excluded because they are in a plate boundary zone)  The artefacts begin to appear as the distance to the tracking stations is approximately equal to or greater than the spatial resolution of the basis function level (Table 22).
 Varying the network distribution does vary the manifestation of artefacts, on the whole these appear in similar areas with a similar magnitude; notable is Alaska which is more susceptible to artefacts at a lower level (possibly due to the reduction of North American stations).
 The full network velocity biased solutions will be used to create a mask highlighting the regions of confidence for the real solutions, section 6.4.

Seasonal Basis Function Estimations
Each week a loading inversion is calculated from which a linear trend is fitted to each coefficient to produce the secular estimation. This however disregards any seasonal variation which in itself is an interesting point of study. This section is a small deviation away from secular loading before returning to the main focus of the study in section (6.4).
Each week a basis function estimation for a particular level is calculated, from this a weekly grid can be produced which is representative of the seasonal  It is clear that in winter there is an excess of mass in the Northern hemisphere which corresponds to the accumulation of surface water, predominantly snow and ice. This trend is reversed in the summer months where there is a mass deficit in the North. This is documented by Lavallee and Blewitt (2002) who find a degree-1 signal. However Figure 59 is based on a very low resolution estimate and in reality there is unlikely to be a uniform mass transfer across the entirety of Europe and Africa; to identify the pockets of seasonal mass variation, higher level must be estimated. Moving to a level two estimation (estimating nine parameters) Figure 59 reveals mass variation are concentrated in smaller areas as the increased level of estimation refines the spatial resolution. As the level of estimation increases, the combination of the non-secular higher degree coefficients paints a picture of what seasonal surface mass exchange is occurring. As the level of estimation is increased the ability to fit the load to finer detail increases. Figure 60 estimate 16 and 25 parameters per epoch respectively. It is apparent that the higher level estimates are beginning to define distinct locations of periodic surface loading on the surface of the Earth. Discrete pockets are appearing around large river basins (Amazon) which are comparable to GRACE results (Crowley et al. (2008)) and tropical regions which receive strong monsoon rainfall (South East Asia). But there are also the same distinct artefacts appearing in similar areas (e.g. Madagascar, East Africa and South East Asia) that were also present in the synthetic data analysis, section (6.2.2).

Secular Surface Mass Loading
By calculating a simple weighted linear regression of each basis function or spherical harmonic coefficient it is possible to isolate the respective secular loading value. It is this information that I will use in this final section to investigate any present day secular loading. I will do this for the case of removing each of the four previously-described global GIA models and for the null GIA model. For each scenario I will evaluate basis functions from level four up to a maximum of level seven, but focus on a level five and six estimation as the synthetic tests have shown that the low level estimations have too broad a spatial resolution and higher level estimations (8-12, not shown) become 163 unstable in several regions. It may be expected to see a residual loading signal around the regions which have a significant modelled GIA velocity, e.g. Hudson Bay, Fennoscandia and Antarctica (this could be real or the result of potential GIA model error). Any other signal could be a result of present day surface loading (secular changes due to thinning/thickening of glaciers/ice sheets and the increased/decreased storage of water in river basins) or by the over-fitting of the basis function coefficients in poorly constrained areas at higher levels causing anomalously high estimates, also highlighted by the synthetic data analysis, section (6.2).
The first case presented is the estimated loading without the introduction of a GIA model, i.e. only plate velocities removed. Following this are the results from applying the four respective GIA models; the aim is to highlight areas of improvement and any present day secular surface mass loading free of GIA and tectonics. 165 Figure 61 shows unloading in high latitude areas which have a history of glaciation (North America, Fennoscandia, Siberia), however as the estimation increases to level 6 there is a similar pattern of artefacts emerging elsewhere to that seen in the synthetic studies in section (6.2). With the inclusion of a GIA model in the velocity model ( Figure 62 to Figure 64), the areas of GIA-related deformation should be accounted for, leaving only present day surface mass loading but with possible inversion artefacts due to poor spatial coverage remaining in certain areas. A confidence mask has been applied; this has been created from the velocity biased synthetic data set covering all areas which showed trends of greater than 40mm/yr for each level considered.   Each GIA model will have a different set of specific site velocities, which in turn will imply a different pattern of secular loading trends. Fitting basis function levels four (Figure 62), level five ( Figure 63) and level six ( Figure 64) should begin to highlight areas which are common artefacts to all models, section (6.2) and what are genuine geophysical signals. In general the artefacts appear in areas which have sparse data coverage; this mostly leaves higher latitude sites as candidates for real geophysical signals.
There is disagreement over continental North America as to the magnitude and spatial distribution of the effects of GIA, also seen in section (5.2.2) and the differences between the two families of GIA models. The residuals to the Peltier model predict large areas of loading over Hudson Bay, for all levels (Figure 62-Figure 64). This is unlikely to be a real signal as it only seems to radiate around Hudson Bay; but is more likely an artefact or a GIA error. I would suggest this is a GIA model error as in the synthetic tests; loading over North America is very well constrained and has shown no tendency to produce artefacts over Hudson Bay. The Schotman models show a mixture of loading and unloading but at a fraction of the magnitude of the Peltier models suggesting a better modelling of the North American GIA effect, presenting an improved picture of the potential present day loading. Greenland and Western Europe appears to be much more robust for all GIA models and levels presented, with Antarctica and Alaska becoming more susceptible to artefacts as the estimated level increases.
Summation of the grid values over selected areas for each model will yield a range of estimated mass change rates which can be compared to other studies.
Assuming that there are real signals in areas which are well constrained, then there is possible secular mass unloading over Alaska/Bering straits, Greenland and the Antarctic Peninsula detectable at different basis function levels. Other methods (Wu et al. (2010)) and other systems such as GRACE (Tapley et al. (2004)) have been used to determine similar trends of surface mass loss. GRACE, as with GPS, cannot distinguish between mass loss due to present day trends and those caused by GIA and therefore GRACE is as dependent upon the choice of GIA model as GPS is.

Predicting Mass Loss from Secular Trends
Through the process described in this thesis it should be possible to detect mass loss in the same regions as those detected by GRACE, through the surface deformation of the Earth.

GRACE Surface Mass Observations
One of the products from the GRACE mission is an estimation of the variation of surface masses. By calculating the trend from the monthly GRACE estimates a picture of secular mass variation begins to emerge. Figure 65 Baur et al (2009) proceeded to correct these values, accounting for GIA by subtracting a GIA model ) from the signal. This produces a GIA-corrected estimate of present day secular mass variation ( Figure 66).   (20) Below, Table 24 to Table 27 shows the estimated mass loss values for the respective regions from fitting a weighted linear regression. The propagation of errors through the inversion process will bias the estimation of secular surface mass loading. To account for this fact an associated error budget for each level has been calculated using the same inversion methods but using the RMS of the velocity biased trends from section (6.2.2) as the input for each level. Although there are no stations located directly in the proposed study area, there is a very strong network in North America to constrain the inversion. However the south coast of Alaska runs very close to an active deformation zone which could contaminate the solution. Alaska has also shown a tendency for higher level estimation artefacts. The size and shape of the Alaskan area makes it more suitable to lower level estimation than the Antarctic Peninsula, section (6.5.2).

Alaska
From the synthetic analysis, Alaska seems to be resilient to gross artefacts up to a maximum level 6 which is reflected in the error budget. Accounting only for tectonic motion there is a small net loss up to a level 6 estimation (-68±31Gt/yr), once a GIA model has been introduced then there are two very distinct    There is a mixture of estimates even at low levels with the introduction of a GIA model. The Peltier estimates of mass loss appear to be ~3 times larger than the Schotman models at level four, but with good agreement within each GIA family. I suggest a potential present day mass loss of the interior, of somewhere in between the two models -24±6Gt/yr, which is significant at 2σ but there is such a scatter of results and as noted in section (5.2.2), Antarctica GIA remains a significant uncertainty. different answer to the other GIA models presented here. A much more detailed discussion and comparison of Greenland is conducted below, section (6.6). 176 main Schotman model slightly more pronounced; these motions are not present in Peltier's model.

Comparisons of Mass Loss Estimates
Discussed above consider the various secular mass estimates over selected area with different GIA models introduced. All have shown potential for observing present day mass loss but some have also shown a susceptibility to inversion artefacts (Alaska and the Antarctic Peninsula). It should be noted that the Antarctic Peninsula is approximately half the size of Alaska and Greenland, and the approach is unable to produce consistent results. To be able to put any confidence into an Antarctic Peninsula mass loss estimate, a much denser station network is required; the POLENET campaign (www.polenet.org) is achieving this but it will take several years before any estimation of secular rates can be made.
The methods used in this study are offered as an alternative to compare with other GPS results (Wu et al. (2010)) or with alternate systems such as GRACE (Tapley et al. (2004) suggesting an increasing in mass loss, with an overall trend of -230±33Gt/yr.
Velicognia also calculates a mass change of -143±73Gt/yr over the entire Antarctic ice sheet. A range of values from a variety of authors spans from -66±49Gt/yr to -248±36Gt/yr for Greenland (Table 28). However the NC1 result agrees with Wu et al. (2010) if we consider the GRACE only values, -68±28Gt/yr and -71±6Gt/yr (Luthcke et al. (2008)). This lower estimate could be explained by the lack of stations within the study area, the choice of GIA model used or, as the region is close to a deformation zone. All of these reasons make it increasingly difficult to estimate consistent results due to inversion artefacts.

Summary
In this chapter I have shown that it is possible to fit a set of modified basis

Chapter 7. Conclusions
The precision of station coordinate estimates from GPS measurements is continually improving. This is due to the on-going efforts to improve processing strategies and models, section (3.3). The computed combination of reprocessed data contains ~80 extra stations (GPS week 1000-1300) and has shown great improvement over the operational data, sections (3.6.6, 3.6.7); this is highlighted by the reduced WRMS, Figure 26. This improvement has allowed the construction of an improved global kinematic solution, section (4.1) which is fundamental to the estimation of plate tectonic rotations, section (2.5.3) and surface mass loading, section (2.3).

Reprocessed Data Series
This study pays careful attention to the quality of the weekly combined solution and the kinematic station estimates, Chapter 3 and Chapter 4. The best measures of quality can be concluded by comparing the reprocessed dataset to its predecessor.  The estimated seven parameter Helmert transformation variations are reduced by ~75% across the board. This reduction is even greater for some AC solutions in the very early weeks (Table 16).
 The most important improvement seen between the operational and reprocessed data sets is the removal of systematic offsets (Figure 27, Figure 28 and Figure 29) and noise ( Figure 30 and Figure 31) which was the primary aims of the IGS reprocessing campaign.
 By only including stations with over 2.5 years of data (following recommendations of Blewitt and Lavallee (2002), with the majority of included station having a much longer data span, station velocities, section (4.4), are now estimated with sub-millimetric formal errors, (Appendix A), even after accounting for an appropriate noise model, section (4.4.2).

Model Estimation
A large section of this work has been the development of an integrated 3D  (2010)), (Figure 43). This allows for the removal of a range of a priori GIA trends from the raw velocities, section (5.2), although these small changes have an important effect upon the final estimates. Shows sensitivity to small regionally correlated velocity errors It should be noted that the GIA models used only represent a 1D Earth model and the introduction of a 3D model could potentially introduce horizontal velocities different by ~1-2mm/yr (Latychev et al. (2005a)), which would influence the tectonic plate estimation, but until a 3D GIA model becomes publically available this influence can only by hypothesised.
The use of a combined GIA and plate tectonic velocity model in place of a linear velocity estimate produces residuals which should be due to present day secular surface mass loading and model error. The final stage of work, Chapter 6, fits a set of modified spherical harmonic functions to the residual at different levels; this method models and recover the remaining trends.
One of the aims of this study was to determine the maximum level of basis function estimable; preliminary investigations showed that higher level inversions (7-12) were prone to spurious signals, so two tests of global inversion quality were carried out. The first used a global  2 dof to test if there were any indication of a fall in inversion quality as the inversion level increased.
Of the levels tested (1-8) there was no indication that the higher levels were indeed a poor fit with all values falling within the expected  2 dof range (Table   23).
The second test used a trend free synthetic data set, section (6.2), to which several scenarios of random white noise (WN) and coloured noise as a velocity bias (VB) were added. The inversion was completed as with the real data sets and the RMS of WN and VB results were produced. The WN test shows that for levels (4-6), with no trend was expected, there is only a very small loading trend of 2-3mm/yr introduced through the inversion process which is deemed to be the random error of the basis function estimation.
For the WN+VB synthesis which has a small velocity bias added (standard deviation of 0.2mm/yr and 1mm/yr for horizontal and vertical components respectively), as the basis function level increases, there is a growing susceptibility to inversion artefacts for regions at different levels. These artefacts tend to appear where the spatial data coverage is sparse, Figure 55, whereas areas which are, in or around, dense GPS network coverage show resilience to inversion artefacts. This synthesis indicates which observed loading signals are potentially real, and provides an estimation of the inversion error. The same VB test but using the reduced network shows similar spatial pattern and magnitude of inversion error, section (6.2.2).
Despite the different GIA models, section (2.4.4), similar areas of present day mass loss become apparent at the maximum reliable inversion level Alaska (level six), Antarctic Peninsula (level four) and Greenland (level seven), Figure   62- Figure 64, with an expected low to zero benchmark test of Australia (level seven). There are some differences between the GIA adjusted estimated load in both spatial distribution and magnitude. Alaska and the Antarctic Peninsula tend to suffer from inversion artefacts at a lower level compared to Greenland which is, in part, due to the latter's size, shape and data coverage. This again reiterates the method's reliance upon global data coverage.
Calculation of the total mass change in the selected areas over all levels produces an initial mixed message of mass loss and gain, but when accounting for the maximum reliable inversion level for the selected areas (Alaska, Antarctic Peninsula and Greenland), the message becomes one of clear mass loss, (Table 24-Table 27), with Australia showing a small mass gain within the inversion error.
These offer an independent assessment of global mass change which can be compared to the NC1 estimation. Greenland has by far the most investigations, and also happens to be the best candidate for this method especially at the higher basis function levels. Apart from the ICE5G-VM2 model, there is good agreement at levels six and seven of a value ~-140±30Gt/yr, with the ICE5G-VM2 model giving an estimated value of -234±30Gt/yr. GRACE, as with GPS, is dependent upon GIA model choice. Table 28 shows the wide ranging values of GRACE mass change estimates over Greenland with some of the ICE5G-VM2 adjusted estimates estimating mass change of -234±24 to -248±36Gt/yr, (Chen et al. (2006c), Velicogna and Wahr (2006a)) which are similar to the NC1 ICE5G-VM2 estimate. Considering the various studies over Greenland which range from -66±49Gt/yr to -248±36Gt/yr the NC1 solutions fit well within this range for three out of the four models used with the fourth (ICE5G-VM2) at the higher end of values.
Comparing values over Alaska at a lower level to other studies we find our estimates are ~50% less, this could be due to the "true" loading being contaminated by inversion artefacts, tectonic deformation or network coverage.
Synthetic tests would suggest that the processes is mainly affected by station coverage, combined with residual GIA as being the cause for this reduced value over Alaska. As for the Antarctic Peninsula the level four values suggests a small mass loss but this quickly becomes biased by artefacts and is highly unstable between levels. There is no real comparable study for the Antarctic Peninsula on its own.

Future Work and Recommendations
The kinematic solution and suite of global secular loading models presented here are a valuable source of data for scientific research and I have only followed one line of investigation. Even with all considerations there are still components of this study which could be further extended. In the future additional GNSS data sets will become available (e.g. GLONASS in some AC solutions). The addition of these data sets will continually increase the length and number of stations available for study, and may reduce technique-specific 183 errors. The improvement of the IGS reprocessing campaign's first round is highlighted by the vast increase in the precision compared to the historic operational processing. Although such a dramatic improvement of any further rounds of reprocessing might not be as apparent. The first recommendation made is that these reprocessing campaigns are a valid and valuable use of resources and that future campaigns should be considered as new model improvements come to light. One of the main limitations of this study is the number of stations available, especially in the southern hemisphere. One of the main improvements is due to the densification of the available network especially in early weeks. It cannot be stressed enough the importance of including previously unused stations which has shown to be a factor throughout this study. This should partly be address through the second round of reprocessing is already in the pipeline and also intends to include: The consideration of vertical station motion in this study has been achieved by employing a limited number of global GIA models which use a 1D Earth model.
The inclusion of these models has been shown to reduce the vertical residuals of the secular loading model. The included models have known deficiencies, especially around Antarctica. The second recommendation is twofold, first to extend the range of GIA models used, focussing upon Antarctica, and second to introduce a GIA model based on a 3D earth model.
This initial kinematic solution assumes a white noise station velocity error model which is insufficient to account for the known correlations between weekly GPS observations, individual tracking stations' monument motion or any other low frequency systematic noise. To account for this a uniform variance scaling factor has been applied to the kinematic solutions uncertainties based on work by Santamaria-Gomez et al. (2011) and Bos et al. (2008). The third recommendation would encourage a more in-depth investigation of individual station coloured noise components using the reprocessed data and a more through implementation of these into the kinematic velocity estimation.
This study has only used data from one spaceborne geodetic technique with a limited set of GIA models, and carried out the inversion using basis functions.
The method has been shown to suffer from inversion artefacts over regions with sparse data coverage; the inclusion of additional sites (e.g. relaxing inclusion criteria) would allow additional stations to aid the inversion. The inclusion of data from other missions such as GRACE or sea-surface altimetry (corrected for steric ocean changes) would increase the redundancy of the solution and strengthen any inferences made about present day surface mass changes. The residual loading harmonic estimation only uses continental basis functions. The inclusion of oceanic basis functions, which account for the redistribution of masses within the ocean due to circulation currents, might improve the estimation of loading slightly. This again would provide a better picture of global mass redistribution.
The methods described in this study are presented as an alternative which only uses deformation data. This is an advantage as it is neither dependent upon a particular ocean model output nor dependent upon assumed covariance between GIA deformation fields. This produces a much clearer and more transparent method, i.e. when using model A the predicted mass change is X Gt/yr with a known ocean effect.