C Mitsakaki

A Mimidou

Simultaneous estimation of secular and episodic crustal motion via geodetic observations

Keywords: GPS Positioning, Deformation Measurements, Dislocations, Corinth Rift

Repeated high quality geodetic observations allow the estimation of the free surface velocity field for a region. Usually, a yearly secular rate is estimated, while the possible episodic motion (seismic slip) due to an earthquake is evaluated via inversion of the geodetic data. The episodic motion influence is subtracted from the total field. The actual region and the respective geodetic observations affected by the seismic event are often assessed by rather vague criteria. This paper deals with an attempt to estimate simultaneously the secular and the episodic two-dimensional crustal motion for a region by means of repeated GPS observations of a geodetic network, carried out for a number of epochs.

Positions of the gulf of Corinth GPS network points: red triangles, blue dots and black dots correspond respectively to the 1st order, 2nd order points (also old Hellenic triangulation network points) and to the ones common with the Central Greece network (another series of GPS campaigns in Greece). Dented segments show the localisation of the main active faults in the area. From [START_REF] Briole | A ten year analysis of deformation in the Corinthian Gulf via GPS and SAR Inteferometry[END_REF]. 3 of [START_REF] Avallone | Analysis of eleven years of deformation measured by GPS in the Corinth Rift Laboratory area[END_REF].

A c c e p t e d M a n u s c r i p t

Table 2

Parameters for the offshore faults of the 1992 and 1995 earthquakes.

The parameters are: Ms the magnitude of the event, Mo the seismic moment, Lat and Long the geodetic coordinates of the fault's center,

Length the length of fault, Strike the fault's azimuth, Dip the fault's dipping angle and Top/Bottom the upper and lower depth of the fault [START_REF] Briole | The Ιtea (Μ ~ 5.9) Earthquake of November 18, 1992. -Characteristics of the Main Shock Inferred from Body Wave and Ground Displacement Analysis[END_REF], [START_REF] Bernard | The M S =6.2, June 15, 1995 Aigion earthquake (Greece): Evidence for low angle normal faulting in the Corinth rift[END_REF], [START_REF] Avallone | Analysis of eleven years of deformation measured by GPS in the Corinth Rift Laboratory area[END_REF] 1992 A c c e p t e d M a n u s c r i p t

Table 3

Secular strain parameters for the two homogenously deforming blocks (North and South). At the same time three fault segments are considered for the episodic motion: a single segment model for the 1992 Galaxidi earthquake and a double one (the upper slab and a lower relaxation segment) for the 1995 Aigion event.

A is the azimuth of maximum strain rate (semiaxis of the strain ellipse) e max .

This work Smoothed Data

North block

Secular rate parameters

All Epochs (1990Epochs ( -2001) ) (3 fault segments for the episodic motion) A c c e p t e d M a n u s c r i p t

The last column is taken from previous work [START_REF] Agatza-Balodimou | Deformation Studies in the Corinthian Gulf via Multi-epoch Analysis of Geodetic Data, 11 th FIG Intern[END_REF], [START_REF] Avallone | Analysis of eleven years of deformation measured by GPS in the Corinth Rift Laboratory area[END_REF]. It represents the secular parameters estimated after the average co-seismic motion for the 1995 Aigion earthquake (considered as a single dislocation segment) is eliminated from the GPS data.

A c c e p t e d M a n u s c r i p t Parameters chosen and slip results for the 1992 and 1995 earthquakes with the a posteriori variance factor o σ ˆof the model.

First case: Only the 1992 event modelled; second case: the 1992 and 1995 earthquakes modelled as single blocks; third case: the 1992 as a single block and the 1995 as a two block dislocation; fourth case: only the 1995 event modelled as a two block dislocation.

*

The first values for the a posteriori o σ ˆof the model refer to the case of using all available GPS data (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001) and the second one when only data from the interval 1990-1997.8 have been used.

**

The shadowed rows refer to results from previous studies [START_REF] Briole | The Ιtea (Μ ~ 5.9) Earthquake of November 18, 1992. -Characteristics of the Main Shock Inferred from Body Wave and Ground Displacement Analysis[END_REF], [START_REF] Bernard | The M S =6.2, June 15, 1995 Aigion earthquake (Greece): Evidence for low angle normal faulting in the Corinth rift[END_REF] for comparison.

A c c e p t e d M a n u s c r i p t

*

The data were acquired within the frame of a project in the vicinity of the Gulfs of Euboea and Corinth (Agatza-Balodimou et al., 2003). 1990 1991 1992 1993 1994 1995.5 1995.8 1997* 2001 First order points Second order points The most pronounced characteristic is the high extension rate over the Aegean Sea.

GPS Campaigns

Thus, the southwestern Aegean moves, relative to Eurasia, toward the SSW at about 30-40 mm/yr [START_REF] Nyst | New constraints on the active tectonic deformation of the Aegean[END_REF]. Normal faulting, organized into subparallel systems distributed over areas of tens or even hundred kilometres wide, is a feature often appearing in actively extending regions on the continents.

One of the most prominent and active features of such a system of faults is the rift along the Gulf of Corinth, with a large seismicity rate. It is an area of high seismic risk, since it is surrounded by several large, densely populated cities, among them Athens. Therefore extensive research work has been carried out for more than a decade.

A c c e p t e d M a n u s c r i p t 4

The Gulf of Corinth is the most rapidly extending rift system in Greece with about 120km length and 30km width and a WNW-ESE trend. It is believed to be active at the present rates since the last 5Myrs (Taymaz et al., 1991), [START_REF] Armijo | Quaternary evolution of the Corinth Rift and its implications for the Late Cenozoic evolution of the Aegean[END_REF]. Most of the active surface normal faults are located on the south side of the gulf (North Peloponnesus) while some south-dipping faults appear on its northern side (Sterea Hellas) [START_REF] Armijo | Quaternary evolution of the Corinth Rift and its implications for the Late Cenozoic evolution of the Aegean[END_REF]. Present day extension rates, estimated from GPS observations, for the eastern part of the gulf are significantly lower (~5 -6mm / yr) than the central-western ones [START_REF] Briole | Active Deformation of the Gulf of Korinthos, Greece: Results From Repeated GPS Surveys Between 1990 and 1995[END_REF], [START_REF] Clarke | Geodetic estimate of seismic hazard in the Gulf of Korinthos[END_REF].

Medium size seismic events of the order of M S ~ 5.5 -6.5 are often the case. 1). Today, beyond these periodically observed points, five permanent GPS stations have been established in the area for continuous monitoring [START_REF] Briole | The Gulf of Corinth CORSEIS" permanent GPS network: framework and perspectives[END_REF].

The key objective of this long-term project is to combine methods and data in order to estimate the seismic hazard in the region. Data both from the periodic reoccupation of the dense geodetic GPS network and the ones derived from the small permanent network may be used for this purpose. More explicitly, not only Often the rapture of the fault occurs at depth and no surface trace of the slip is visible. Thus, in order to evaluate episodic motion parameters using the surface velocity field derived from geodetic observations (today mostly GPS ones) one has to solve an inverse problem which has no unique solution. For medium size seismic events a typical procedure is a trial-and-error one, while observations and data from other disciplines (seismology, geology etc) provide constraints for the solution of the problem via a dislocation model (e. g., [START_REF] Snay | Regional deformation of the earth model for the San Diego region, California[END_REF][START_REF] Hudnut | Co-Seismic Displacements of the 1994 Northridge, California, Earthquake[END_REF].

Thus, geodetic network points with repeated observations (e.g., GPS) considered relatively close to the location of the earthquake and therefore probably affected by the co-seismic displacements are chosen. The coordinate time series of these are computed and examined in order to assess whether their positions are influenced from the event and to what extent. A dislocation model is chosen, with some of its parameters constrained by other than geodetic data (e.g., aftershock distribution at depth) and its displacement field on the earth's free surface is compared with the one available from the geodetic observations. The procedure is repeated until the estimated dislocation model parameters best fit the surface displacement field.

In the present work a two-dimensional model is described that solves for both the secular and some parameters of the episodic motion and accepts all available geodetic data (coordinate time series) as observations. Medium to large seismic and Aigion 1995 earthquakes and the results are discussed.

The model

Until today, the common technique for deformation studies using geodetic observations is to carry out campaigns where discrete network points on the earth's surface are occupied repeatedly. Therefore, the data are sparse both spatially and temporally, and they should be homogenized by referring them to the same geodetic reference frame (e.g., ITRF2000) and epoch in order to be further used.

Moreover, geodetic observations, like all types of observations have random errors, while they may be, also, contaminated by systematic errors (biases) 1 or gross errors (mistakes or outliers) 2 . Therefore, their reliability has to be statistically tested. The deformation model that is estimated by the geodetic observations has also to be tested for its congruency and its parameters for their statistical significance. Only the parameters found statistically significant should remain in the final expression of the model.

1 Systematic errors may appear in observations when necessary reductions are omitted, or the influence from various sources (such as the ionosphere effects on GPS measurements) may be not properly modelled etc. They usually obey physical laws and in most cases may be modelled by mathematical expressions.

2 Outliers are mistakes occurring by negligence. They are large in size and relatively easily detected.

A c c e p t e d M a n u s c r i p t 7 So far, the geodetic GPS observations available in Greece, remain mostly sparse and thus, inadequate to discern a combination of parameters regarding the fault geometry. Thus, the detection of which fault is responsible for a seismic event, its geometry and location are information derived primarily from geological and seismotectonic observations. At present, the model described here solves only for the strike and dip slip components of the episodic motion of a seismic event that occurred during the time interval for which adequate GPS observations are available. The software uses as data the displacement at a given position if that location is within the specified radius of the centroid of the earthquake.

With respect to the secular motion the area is considered as consisting of a mosaic of blocks with borders active fault traces or lines that distinguish changes in the velocity field.

Since the accuracy of the vertical component determined from GPS observations is at least two to three times worse than the horizontal components the twodimensional approach was chosen in this study. Thus, a two-dimensional model was developed that solves for both the secular and the two episodic motion parameters (strike and dip slip). Due to the average size of displacements, the infinitesimal elastic strain theory is followed for the secular motion without loss of the required precision. The data consists of the two-dimensional projection coordinates of the network point positions (E, N) for the various epochs of observations derived from the three-dimensional GPS (X, Y, Z) ones. Often, maximum and minimum strain rates (i.e. the axes of the strain rate ellipse) e max and e min are also used, as well as the total shear rate used in the model which is defined by the conditions [START_REF] Snay | Regional deformation of the earth model for the San Diego region, California[END_REF]: Conditions (2) consider that the slip on the m-th parallelogram occurs instantaneously at the instant t m even though this may result in mistakenly assigning a post-seismic activity to the co-seismic phase. However, a more refined model that would attempt to distinguish between the phases of the seismic cycle (pro-, co-and post-seismic cycle) would require continuous monitoring observations, at several network points. This, until recently, was not easily realized in the case of the GPS networks and still remains a demanding and rather expensive approach. defining the origin of a local three-dimensional right-handed Cartesian system (x 1 , x 2 , x 3 ). The axis x 1 coincides with the direction of the fault dip, while the axis x 2 is taken parallel with the fault strike (Figure 3).

) γ + γ ( = γ 2 2 2 1 1/2    .                                                     m m m 0 0 o m m a m j 0 j 0 j i 0 j i k 1 2 2 1 0 j k N E k 0 i j i 0 i j i ds ss ) N , E ( F - N) (E, F R ) t , r(t ) t - t ( N - N E - E 2 γ - ρ ) 2 γ - ω ( - 2 γ + ω 2 γ + ρ + ) t - t ( d d = N - N E - E           A c
for t m < t o         m j m j m j t t 0 t t 1 - ) t , r(t (2) for t m > t o         m j m j m j t t 1 t t 0 ) t , r(t
The displacement, in map projection coordinates, of an arbitrary point

) N , E ( j i j i
on the earth's surface due to the m-th parallelogram is analysed in two components.

One is oriented parallel to the fault azimuth (along strike)

N) (E, u 1m

and the other transversely to this

N) (E, u 2m

. Then, from dislocation theory there are functions f im (E, N) and g im (E, N) (with i =1, 2) such that:

(3)

Functions f im and g im depend on the point's position. They, also, depend on the geometry, the position and the orientation of the parallelogram as well as the Poisson ratio ν, which characterizes the elasticity of the earth's crust. In this model the Okada dislocation formulae are used [START_REF] Okada | Surface deformation due to shear and tensile faults in a half-space[END_REF].

Formulae (3) refer to the local three-dimensional system (x 1 , x 2 , x 3 ), with origin the centre of the fault trace at the free surface. The origin's coordinates in the map

                    m m 2m 2m 1m 1m 2m 1m ds ss N) (E, g N) (E, f

N) (E, g N) (E, f N) (E, u N) (E, u

A c c e p t e d M a n u s c r i p t 12 projection reference system are E C , N C . This point is considered free of displacement.

For a point on the earth's surface, with map projection coordinates (E_, N_) before the m-th earthquake, the relation that connects the two-dimensional geodetic reference system (E, N) with the arbitrary local reference system (x 1 , x 2 ) is:

(4)
where α is the fault azimuth (strike).

The contribution of the dislocation expressed in the local system (x 1 , x 2 ) due to the m-th event is the vector

) u , (u 2m 1m
. Then, directly after the earthquake, the new projection coordinates of the point (E + , N + ) are expressed as:

(5) where the point (E C , N C ) is considered as free from motion.

Combining expressions (4) and ( 5) the influence of the m-th dislocation on the map projection coordinates of the point is:

(5')

The motion due to the earthquake at the network reference point according to (5'). When the motion of the arbitrary point (Ε i , Finally, the formula (1) for each k -th block and all m seismic events may be generalized as:

                     2 1 C C x x cosα sinα sinα cosα N - Ν_ E - Ε_                          2m 2 1m 1 C C u x u x osα sinα sinα cosα Ν - N Ε - E c                               
(7)
The variables in formula (7) are:

                 m m 0 0 o m m a m m j ds ss ) N , E ( F - N) (E, F R ) t , r(t t δ ) N , E ( δ E + d = δ 0 0 k k                                o 2m 2m o 1m 1m a m m u u u u R N Δ E Δ o _ o o _ o Ν - N Ε - E                            m m o o o m m a o _ o o _ o m m ds ss ) N , E ( F - N) (E, F R Ν - N Ε - E N Δ E Δ        N) (E, g N) (E, f N) (E, g N) (E, f N) (E, F 2m 2m 1m 1m m        ) N , E ( g ) N , E ( f ) N , E ( g ) N , E ( f ) N , E ( F o o 2m o o 2m o o 1m o o 1m o o o m
A c c e p t e d M a n u s c r i p t

14 k 0 i j i 0 i j i N - N E - E =         δ , k d  = k N E d d           , k E  = k 1 2 2 1 2 γ - ρ ) 2 γ - ω ( - 2 γ + ω 2 γ + ρ                       , ) N , E 0 0 ( δ =         N - N E - E 0 j i 0 j i
The method described here has similarities to block modelling (e.g., [START_REF] Mccaffrey | Crustal Block Rotations and Plate Coupling, in Plate Boundary Zones[END_REF], although ours is a much simpler model. For the part of the surface velocities due to a locked fault (episodic motion) [START_REF] Mccaffrey | Crustal Block Rotations and Plate Coupling, in Plate Boundary Zones[END_REF] uses the same dislocation formulae [START_REF] Okada | Surface deformation due to shear and tensile faults in a half-space[END_REF]. However, the approach in McCaffrey ( 2002) refers to block-bounding faults that follow the stick-slip model (aseismic steady state slip and slip 'deficit', as it is called in the locked phase). The present work is better suited for active faults, not of necessity bordering a block, that slip due to an earthquake.

In the present work, the secular motion is described by the infinitesimal strain rate inside each block. We preferred strain rate modelling, instead of the GPS velocity field, since the requirement for a strictly common spatio-temporal reference frame is not essential in our case. Part of the strain rate tensor is the solid body rotation rate of the block. This rotation rate may be transformed to an Euler rotation rate according to [START_REF] Mccaffrey | Crustal Block Rotations and Plate Coupling, in Plate Boundary Zones[END_REF].

Another difference lies in the data used in the present work. The data are positional information from several GPS campaigns, and the model solves for strain A c c e p t e d M a n u s c r i p t 15 rates and co-seismic slips (formula 7). The data in the block modelling approach is, usually, smoothed averaged velocities that rather lack the inconsistencies of real data. However, they still need to refer to the same geodetic reference frame and the same epoch (e.g., ITRF2005, epoch 2007) in order to be used as data for deformation analyses [START_REF] Mccaffrey | Crustal Block Rotations and Plate Coupling, in Plate Boundary Zones[END_REF], [START_REF] Nyst | New constraints on the active tectonic deformation of the Aegean[END_REF].

Analysis

Since 1990 the Higher Geodesy Laboratory of the National Technical University of Athens participated in a European multi-disciplinary research programme for monitoring the tectonic behaviour of the Corinth rift region. A GPS network, part of a much larger network was established in Greece by several research teams [START_REF] Briole | Active Deformation of the Gulf of Korinthos, Greece: Results From Repeated GPS Surveys Between 1990 and 1995[END_REF]. The network includes about 50 first order points (measured at least three times in a given campaign). In addition, about 150 second order points were observed one or two times during at least one campaign. It should be mentioned that approximately 140 of these GPS network points are pillars of the Hellenic triangulation network.

Eleven GPS campaigns were carried out from 1990 to 2001 (Table 1). Two of them (November 1992 and June 1995) took place after the M S = 5.9, 18 November 1992, Galaxidi and the M S = 6.2, 15 June 1995, Aigion earthquakes.

A c c e p t e d M a n u s c r i p t 16 For all campaigns, IGS precise orbits and data from IGS stations were used to tie the network to ITRF2000. No ambiguity fixing was allowed for baselines longer than 100 km, in other words between points in the network and the IGS sites. The average percentage of ambiguities fixed in the network was ~ 85% [START_REF] Avallone | Analysis of eleven years of deformation measured by GPS in the Corinth Rift Laboratory area[END_REF]. Expected sources of errors in each campaign's final results are the centring of the antennas for the horizontal components, and the antenna heights and troposphere modelling for the vertical component [START_REF] Avallone | Analysis of eleven years of deformation measured by GPS in the Corinth Rift Laboratory area[END_REF]. So far, only horizontal solutions of the GPS data have been employed for the Gulf analyses; the same applies here. Time series of the map projection coordinates offer average uncertainty estimation of the order of a few mmyr -1 , which is a more realistic estimate than the uncertainties based on the GAMIT solutions [START_REF] Avallone | Analysis of eleven years of deformation measured by GPS in the Corinth Rift Laboratory area[END_REF].

Previous extensive research work carried out for the Gulf of Corinth shows that the extension of the rift is localised along a narrow offshore zone of about 10km

width. The extension rate, as derived from about 11 years of GPS observations is not uniform, with an average of 11mm/yr in the central part of the rift (largest rate of 16mm/yr close to Aigion town) and tapering off to 5-6mm/yr at the eastern edge of the gulf [START_REF] Avallone | Analysis of eleven years of deformation measured by GPS in the Corinth Rift Laboratory area[END_REF], [START_REF] Agatza-Balodimou | Deformation Studies in the Corinthian Gulf via Multi-epoch Analysis of Geodetic Data, 11 th FIG Intern[END_REF], [START_REF] Briole | A ten year analysis of deformation in the Corinthian Gulf via GPS and SAR Inteferometry[END_REF].

A c c e p t e d M a n u s c r i p t 17 In view of these findings and in order to keep the model as simple as possible it was decided to consider the area of the gulf as consisting of two blocks, the north and south.

It is estimated that for the central part of the rift all recent large earthquakes (Eratine of Phokida, M S = 6.3, 1965; Antikyra, M S = 6.2, 1970; Galaxidi, M S = 5.8, 1992, Aigion, M S = 6.2, 1995) activated offshore faults with shallow north-dipping planes [START_REF] Baker | Earthquake mechanisms of the Adriatic Sea and Western Greece[END_REF]. In contrast the deformation pattern of the western part of the Corinth rift differs in that all active normal faults dip at large angles (50° to 60°) suggesting a probably significant structural change [START_REF] Bernard | Seismicity, deformation and seismic hazard in the western rift of Corinth: New insights from the Corinth Rift Laboratory (CRL)[END_REF].

Since the first GPS campaign took place in 1990, it was decided, at first, to investigate only the Aigion 1995 earthquake for the episodic motion. The event was large enough to affect the displacement field of a broad region and GPS data were available both before and after its occurrence.

The Galaxidi 1992 earthquake, an event of smaller size, took place in 1992. At the time, the GPS network was only partially established. Thus, the velocity field was poorly resolved. However, an attempt was made to include this event in the analysis.

The parameters for the two earthquakes taken from previous studies are presented in Table 2 [START_REF] Bernard | The M S =6.2, June 15, 1995 Aigion earthquake (Greece): Evidence for low angle normal faulting in the Corinth rift[END_REF], [START_REF] Briole | The Ιtea (Μ ~ 5.9) Earthquake of November 18, 1992. -Characteristics of the Main Shock Inferred from Body Wave and Ground Displacement Analysis[END_REF], [START_REF] Mitsakaki | Coulomb Stress Changes in the Gulf of Corinth (Greece) for the 1992-1995 Period[END_REF].

Secular motion parameters

The most characteristic cases for the secular motion parameters of the North and South blocks are depicted in Table 3. It should be mentioned here that the displacement residuals were statistically inspected for outliers at the 99% confidence level. The respective suspect observations were eliminated and the revised data set was used for a new solution. The last column (Table 3) is extracted from previous work [START_REF] Agatza-Balodimou | Deformation Studies in the Corinthian Gulf via Multi-epoch Analysis of Geodetic Data, 11 th FIG Intern[END_REF], [START_REF] Avallone | Analysis of eleven years of deformation measured by GPS in the Corinth Rift Laboratory area[END_REF] and presented here for comparison. The effect of the average co-seismic motion for the 1995 event was estimated and eliminated from the velocity field [START_REF] Avallone | Analysis of eleven years of deformation measured by GPS in the Corinth Rift Laboratory area[END_REF]. Then the secular parameters were estimated [START_REF] Agatza-Balodimou | Deformation Studies in the Corinthian Gulf via Multi-epoch Analysis of Geodetic Data, 11 th FIG Intern[END_REF].

In the present work the secular parameters of the 1990-1997.8 dataset appear to be in better agreement with the smoothed data ones -at least for the north block.

This may be due to the way the smoothed field was estimated by Avallone et al. (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001). A co-seismic offset for 15 points, located close to the 1995 event's epicentre, was estimated from the regression models. This co-seismic effect was subtracted from the total field of the affected points. Then the velocities estimated from the regression models were considered as the smoothed secular velocity field.

Since the north block GPS points were the ones mostly affected by the 1995 event,

the removal of the co-seismic effect by [START_REF] Avallone | Analysis of eleven years of deformation measured by GPS in the Corinth Rift Laboratory area[END_REF] is probably the critical factor for the smoothing out of the remaining secular displacement field.

This smoothed field was used in [START_REF] Agatza-Balodimou | Deformation Studies in the Corinthian Gulf via Multi-epoch Analysis of Geodetic Data, 11 th FIG Intern[END_REF] but the two blocks were treated separately and the respective secular parameters were independently estimated.

In the present work not only the secular motion for the two blocks but also the episodic motion was estimated simultaneously. Therefore, the effect of the 1995 event's episodic motion was also relieved in our model, especially for the north block. Hence, the secular motion solution of [START_REF] Avallone | Analysis of eleven years of deformation measured by GPS in the Corinth Rift Laboratory area[END_REF] of the smoothed field for the north block is very similar to ours.

In contrast, the present solution for the south block is similar whether the 1990-1997.8 dataset or all data (1990 -2001) are used (Table 3, 1 st and 2 nd column).

These secular strain rates are quite lower compared to the ones estimated in the previous work (Table 3, last column) [START_REF] Agatza-Balodimou | Deformation Studies in the Corinthian Gulf via Multi-epoch Analysis of Geodetic Data, 11 th FIG Intern[END_REF][START_REF] Avallone | Analysis of eleven years of deformation measured by GPS in the Corinth Rift Laboratory area[END_REF]). It appears that the present model allows for some network points of the south block (located around the region of Aigion and close to the coast) to be partially influenced by the co-seismic effect of the 1995 earthquake. Therefore, a

Episodic motion parameters

In all cases the geodetically derived co-seismic fault slip was significantly larger than the values derived in previous studies (Table 4). Several trials took place in order to estimate the model that better fitted the data. Global congruency testing for choosing the best fitting model for the Galaxidi and Aigion events was carried out for all cases. No weighting of the GPS data was considered; therefore the test statistic used was: For the 1992 Galaxidi event, GPS data from the period 1990-1994 were used.

The estimated slip is quite larger (~120cm total slip) than the previously derived one of 12cm (Table 4, first case, 2 nd row) [START_REF] Briole | The Ιtea (Μ ~ 5.9) Earthquake of November 18, 1992. -Characteristics of the Main Shock Inferred from Body Wave and Ground Displacement Analysis[END_REF], while the two slip 1985). Thus, high rates of deformation parameters may still be statistically insignificant if properly weighted.

Furthermore, [START_REF] Briole | The Ιtea (Μ ~ 5.9) Earthquake of November 18, 1992. -Characteristics of the Main Shock Inferred from Body Wave and Ground Displacement Analysis[END_REF] used not only the relatively few GPS data available at the time (1990)(1991)(1992) but also seismotectonic information to constrain some of the dislocation parameters. The present study used only data from the GPS campaigns. Thus, their dislocation positional parameters (Table 1) were chosen as known and only the strike and dip slips were estimated. It appears, though, that the geodetic data alone may not be sufficient to constrain the episodic motion for the Galaxidi event properly.

The slip for the 1995 earthquake (Table 4, 8 th row) is also larger than the 87cm estimated by [START_REF] Bernard | The M S =6.2, June 15, 1995 Aigion earthquake (Greece): Evidence for low angle normal faulting in the Corinth rift[END_REF]. However, a slip of this size is in accordance with the Harvard solution of Mo = 5.1x10 18 Nm [START_REF] Bernard | The M S =6.2, June 15, 1995 Aigion earthquake (Greece): Evidence for low angle normal faulting in the Corinth rift[END_REF].

Alternatively, using the seismic moment (Table 2, 2 nd row) given by [START_REF] Bernard | The M S =6.2, June 15, 1995 Aigion earthquake (Greece): Evidence for low angle normal faulting in the Corinth rift[END_REF], and a slightly lower rigidity of μ = 2.9x10 10 Nm -2 , instead of μ = 3.39x10 10 Nm -2 used in their study, the slip is again of the order of 1m. Besides, the For a number of network points, the size of the remaining residuals, after the 1992 and 1995 earthquakes had been modelled, was of the order of several cm (Table 4, second case). These points were rather far away from the 1995 Aigion earthquake area. This behaviour indicated possible post-seismic relaxation for the 1995 event. In order to evaluate whether such an event took place, the 1995 Aigion earthquake was modelled as a two block dislocation model. The upper one describes the co-seismic episodic motion, as previously discussed. The deeper one, with the same position and strike angle as the upper one, deals with the postseismic relaxation (Table 4, third and fourth case, 3 rd and 2 nd rows respectively) [START_REF] Ergintav | Post-seismic Deformation near the zmit Earthquake (17 August 1999, M 7.5) Rupture Zone[END_REF], [START_REF] Kurt | Geodetic observations of postseismic transients in the context of the earthquake deformation cycle[END_REF].

Global congruency testing indicated that the best fitting model for the Aigion event was the one with two blocks. The statistical significance of the model parameters was also tested. In all two-block trials the model parameters for the episodic motion were found significant at the 95% confidence level.

Finally, the two-block model for the Aigion 1995 earthquake (i.e., Table 4, third and fourth case, 3 rd and 2 nd rows respectively) was chosen as best fitting the data.

The upper block has the geometry estimated from [START_REF] Bernard | The M S =6.2, June 15, 1995 Aigion earthquake (Greece): Evidence for low angle normal faulting in the Corinth rift[END_REF] A statistical value that merits some reflection is the a posteriori o σ ˆof the deformation model (Tables 3 and4). The inclusion of more epochs results in increasing the number of observations which, in turn, raises the degrees of freedom.

Conclusions

The model described in this paper appears to be able to estimate the parameters for both the secular and episodic motion for a number of blocks and seismic events.

However, a displacement field with good temporal and spatial coverage is necessary.

A modification of the model to account for possible post-seismic relaxation may improve its sensitivity to better discern between the inter-seismic and post-seismic behaviour of the area under study. The implementation of any post-seismic relaxation model will, obviously, increase the number of unknown parameters. The complexity of such a model would need data from continuous GPS stations in order to better resolve the overall tectonic behaviour of an area.

Furthermore, a full three-dimensional model using the GPS coordinates (X, Y, Z) should be considered in order to describe more accurately the actual threedimensional dislocation model.

Finally, a topic that merits investigation is whether full covariance matrices for the coordinates would provide better estimates for the reliability of the model and the statistical significance of its parameters.

  Figure 1Positions of the gulf of Corinth GPS network points: red triangles,

Figure 2

 2 Figure 1

  Figure 2

  the inter-seismic strain rate across the Gulf of Corinth but also the estimation of the geometric parameters (co-seismic fault slip) of the episodic motion associated with seismic events is sought after.

  velocity field of broad areas. The radius of interference is represented by a constraint embedded in the computer programme and easily adapted to the size of the earthquake. The model is applied for the Galaxidi 1992

  written and used in previous studies (e.g.[START_REF] Agatza-Balodimou | Deformation Studies in the Corinthian Gulf via Multi-epoch Analysis of Geodetic Data, 11 th FIG Intern[END_REF], was modified in order to solve for the episodic motion parameters with their respective full covariance matrices so as to test their statistical significance. It should be mentioned that the displacement residuals may be statistically inspected for outliers at a certain confidence level. Suspect observations are signalled by the programme but their elimination depends on the software user. The software accepts data consisting of as many files as the available epochs of observations. Each file referring to epoch t j contains the arrays of the projected plane coordinates derived at the respective epoch. Furthermore, the program allows for breaking up the region into several blocks. The solution of the parameters is carried out using simultaneously all the pertinent data files for each block and/or seismic event. The advantage of the program lies in its ability to utilize positional information sparse both in space and time and solve for the deformation rate inside each block as well as the episodic motion.The map projection coordinates (Ε i , Ν i ) of a network point i, belonging to the kth block, at observation epoch t j , are modelled as function of the following; its coordinates (Ε ο , Ν ο ) at reference epoch t o , the secular homogeneous infinitesimal strain rate and the contribution of the m-th seismic event that occurred at t m time[START_REF] Mitsakaki | A Method for Deformation Studies in Two Dimensions Using Geodetic Data[END_REF]: body translation and the secular homogeneous infinitesimal strain rates (first two terms of formula 1) the variables are: Easting and Northing of point i at observation epoch j Easting and Northing of point i at reference epoch t Easting and Northing of the network reference point, considered stable both in space and time

  c e p t e d M a n u s c r i p t 10 The last term in formula (1), expressed as a summation, refers to the episodic motion. The term describes the change in position ) the i-th point due to the slip components of the slipping fault (ss m and ds m : slips on the fault surface parallel and transverse to the fault strike respectively). This change takes place at the instant t m of an earthquake occurring on the m-th parallelogram that represents the fault surface. The total effect on the position of each point is given by summing up all contributions from the m-parallelograms.The dependence of the episodic motion on time is made obvious in the case of the step function )

  the last term of expression (1) correlate the fault slips ss m and ds m with the displacements on the earth's free surface. Quite often the fault rupture due to an earthquake does not reach the surface. If the fault rectangle is buried, the upper edge of the parallelogram that simulates the fault surface is projected vertically until it meets the earth's surface. The middle of this line is regarded as the centre of the fault )

  network reference point.

  carried out with respect to the sensitivity of the model to various dislocation parameters. Thus, different depths and dip angles were used and the resulting strike and dip slips were compared. Usually, the model appeared to be insensitive to minor changes of the aforementioned parameters. Radical changes of depth and dip angles affected the model but also decreased the statistical reliability of the solution.

  . They used linear regression models making use of all the GPS time series

  displacement field for this block remains, providing smaller secular strain rates.

σ

  In the formula υ is the vector of the model residuals, n the observations and m the parameters to be estimated), of the Fischer distribution. The statistical significance of the model parameters was also tested using the test statistic ˆis the standard error of the i x ˆparameter and r is the degrees of freedom of the model.

  significant at the 95% confidence level (Table4, first case, 1 st row). However, these high rates may not represent realistically the true motion. This discrepancy may be due to the fact that no full covariance matrices were used for weighting the geodetic coordinates in the deformation model. The use of the full variance -covariance matrix of the coordinates as the weight matrix for a least squares adjustment of a deformation model controls the statistical significance of all the parameters. This significance has to do not with their high value but their real contribution to the model (Agatza-Balodimou and Μitsakaki,

  . (1997) solution was derived using not only the GPS data but also the InSAR images available for the region.

  , with a total slip of 100  12cm. The deeper block dips 80 degrees, extends from 9.5km to 20km depth and has a length of 20km. In this case the a posteriori o σ ˆof the deformation A c c e p t e d M a n u s c r i p t 23 model improved significantly (i.e., from 0.038m to 0.022m in the case of using all epochs).When data from the 1990-1997.8 interval were used the a posteriori o σ was improved from 0.030m to 0.012m respectively.

  However, the aposteriori o σ ˆof the model shows higher values in this case. The network around the Gulf of Corinth was established in a piecewise manner and the network geometry (points observed) differed in each campaign. Therefore, the remaining undetected errors (systematic errors and/or outliers) affect differently the precision of the respective coordinates of each campaign. In other words, not all epochs of observations are of equivalent quality. Since the model used these coordinates as data the residuals of the model are influenced. More explicitly, using more epochs may increase the degrees of freedom but, also, adds to the size of the quadratic sum of the residuals. Hence, the larger values of the aposteriori o σ ˆfor the model. A last comment on the values of the aposteriori variance of the model; the aposteriori o σ ˆof the two-block model for the Aigion earthquake is clearly improved against the single block model. This is evident even in the case discussed previously, when data from all epochs were used (Table 4, last column, first values of o σ ˆ).

and 1995 Earthquakes

	Date	Ms	Mo	Lat	Long	Length	Strike	Dip	Top/Bottom
			(Nm)	degrees	degrees	(km)	degrees	degrees	(km)
	18/11/1992 5.9 0.5×10 18 38.30	22.45	14	270	30	5.2/9.7
	15/06/1995 6.2 3.9×10 18 38.36	22.20	15	277	35	4.5/9.7

Table 2

 2

Epochs 1990-1997.8

				(3 fault segments for the
				episodic motion)
		E d  (mm/yr)	-2±1	-3±2	0±1
		N d  (mm/yr)	5±2	8±2	8±1
		ρ  (μstr/yr)	.12±.02	.11±.03	.11±.03
		ω  (μstr/yr)	.08±.02	.11±.03	.10±.03
		γ (μstr/yr)	.30±.02	.22±.03	.22±.03
		Α (degrees)	17 o .2	17 o .1	13 o .6
		max e  (μstr/yr)	.25±.03	.22±.04	.21±.04
		min e  (μstr/yr)	-.02±.02	-.00±.02	-.00±.02
		E d  (mm/yr)	-0±1	-1±1	-2±1
		N d  (mm/yr)	1±1	2±1	-1±1
		ρ  (μstr/yr)	.02±.02	-.00±.02	.08±.02
	South block	ω  (μstr/yr) γ (μstr/yr) Α (degrees)	.02±.02 .14±.02 16 o .1	.03±.02 .10±.03 16 o .9	.01±.02 .21±.02 5 o .3
		max e  (μstr/yr)	.10±.03	.04±.03	.18±.03
		min e  (μstr/yr)	-.05±.02	-.05±.02	-.03±.02
		o σ ˆ (mm/yr)	±22	±12	±3

Table 3

 3

Slip for the 1992 and 1995 Earthquakes

	Cases examined				Length (km)	Dip degrees	Top/Bottom (km)	Strike slip (m)	Dip slip (m)	o ˆ* σ (m)
	First case	this study Slip values from (Briole et al.,	Only 1992 Galaxidi	14 14	30 30	5.2/9.7 5.2/9.7	.95 ± .06 .90 ± .08 Total slip 0.12m	±.02
		1993)**
	Second		1992 Galaxidi	14	30	5.2/9.7	.92 ± .02 .93 ± .05
	case	this study	1995 Aigion (single block)	15	35	4.5/9.7	1.0 ±.03 1.0 ± .04	±.038 / ±.030
	Third		1992 Galaxidi	14	30	5.2/9.7	.95 ±.01 .91 ± .03
	case	this study	1995 Aigion (2 blocks)	Upper block Lower block	15 20	35 80	4.5/9.5 1.01 ± .02 1.0 ± .02 9.5/19.5 1.02 ±.01 .70 ±.025	±.022 /±.012
	Fourth case	this study	1995 Aigion (2 blocks)	Upper block Lower block	15 20	35 80	4.5/9.5 9.5/19.5 1.04 ± .02 .68 ± .03 1.0 ± .02 1.0 ± .02	±.022 / ±.012
		Slip values
		from (Bernard et al., 1997)	Only 1995 Aigion (Single block)	15	35	4.5/9.7	Total slip 0.87m
		**
	1

Table 4 A

 4

c c e p t e d M a n u s c r i p t

Table 4

 4

Table 1

 1 Overview of the GPS campaigns for the Gulf of Corinth.

Acknowledgments

The authors would like to thank Prof. M. Sakellariou for offering valuable advice.

This work was supported by the Basic Research Programme Leykippos (2006Leykippos ( -2008) ) of the National Technical University of Athens.


Sabine Rödelsperger
email: roedelsperger@geod.tu-darmstadt.de

Matthias Becker

Carl Gerstenecker

Gwendolyn Läufer

Kathrin Schilling

Dieter Steineck

Digital Elevation Model with the Ground-based SAR IBIS-L as Basis for Volcanic Deformation Monitoring

Keywords: Ground-based Synthetic Aperture Radar (GB-SAR), IBIS-L, Digital Elevation Model (DEM), Atmospheric effect, Volcano Fast Response System

Within the project Exupéry, a hybrid deformation observation system is developed, which is part of a Volcano Fast Response System (VFRS). The VFRS shall be a mobile, robust, realtime, easy to use early warning system, which, in case of a volcanic crisis, will support the locale authorities in their decisions about hazard mitigation provisions, especially about the evacuation of people.

The hybrid deformation observation system combines a ground-based Synthetic Aperture Radar named IBIS-L with a network of GPS receivers and will allow the continuous, weather independent determination of areal 3D displacements.

In this paper we present the results of first tests with IBIS-L in an active quarry near Dieburg, Germany. A Digital Elevation Model (DEM) was determined with IBIS-L and compared with a DEM derived by terrestrial laserscanning and photogrammetry. The ability to determine accurate DEMs with IBIS-L is the basis for the combination of different observation techniques, as without a DEM the displacements observed with IBIS-L cannot be georeferenced.

The standard deviation of the differences between the DEM by IBIS-L and the DEM by laserscanning and photogrammetry depends on the surface characteristics. We found a standard deviation of 0.8 m at a slope of rocks and debris, 2.0 m at mining terraces with low vegetation and 3.0 m in a vegetation covered area.

Introduction

Within the framework of the German research program "Geotechnologien", funded by the Bundesministerium für Bildung und Forschung (BMBF), the multi-disciplinary project Exupéry has the aim to develop a mobile Volcano Fast Response System (VFRS) for hazard mitigation (see http://www.geotechnologien.de and http://www.exupery-vfrs.de).

The system allows a quick deployment at any volcano in case of indications for a volcanic crisis or volcanic unrest. The core of the VFRS consists of well-known volcanic monitoring techniques, both space and ground based, to observe seismicity, ground deformation, gas emissions and thermal activity. All data are processed in near real-time and stored in a central data base. Models will be derived to determine the actual state of activity and to set alert levels. The different tasks of Exupéry are described in more detail in [START_REF] Zaksek | Exupery: The mobile Volcano Fast Response System[END_REF].

Within the Exupéry project, the Institute of Physical Geodesy, Technische Universität Darmstadt, Germany develops a hybrid ground based sensor system to monitor displacements of the volcano surface. Such a system has to fulfil certain requirements as:

-high mobility to enable a quick deployment of the system -weather independent continuous and automated monitoring -high temporal and spatial resolution -real-time evaluation Classical approaches for monitoring ground deformation (e.g. extensometers, inclinometers, tiltmeters) require the installation in the volcanic area which might be inaccessible during a volcanic crisis. Therefore the system consists of a network of dual and single frequency GPS receivers in combination with a Ku-band ground-based Synthetic Aperture Radar (GB-SAR) named IBIS-L. The low-cost single-frequency GPS receivers will be stationed inside the radar beam, while the precise dual frequency receivers form the base network. In a new processing approach a joined analysis of discrete GPS data and areal interferometric SAR data, will Many studies show the potential of GB-SAR in applications like displacement monitoring at landslides (e.g. [START_REF] Luzi | Advances in ground-based microwave interferometry for landslide survey: a case study[END_REF][START_REF] Antonello | Ground-based SAR interferometry for monitoring mass movements[END_REF], avalanches (e.g. [START_REF] Martinez-Vazquez | Monitoring structural changes and Stability of the snow cover with a Ground-based Synthetic Aperture Radar[END_REF] and glaciers (e.g. [START_REF] Noferini | Monitoring of Belvedere Glacier using a wide angle GB-SAR interferometer[END_REF]. The advantage compared to optical remote sensing instruments is the ability to monitor areal displacements at any weather condition independent on daylight with a measurement accuracy of 0.1 to 1 mm.

A further application of GB-SAR is the generation of a Digital Elevation Model (DEM) (e.g. [START_REF] Pieraccini | Terrain mapping by ground-based interferometric radar[END_REF][START_REF] Noferini | DEM by Ground-Based SAR Interferometry[END_REF]. The possibility to generate a DEM with IBIS-L is very important within the Exupéry project. The goal of the project is to combine several observation techniques to determine the actual condition of a volcano, which necessitates a common reference system. A DEM is required to georeference the measured displacements with IBIS-L, i.e. to project the displacements into the common reference system. As result of the Shuttle Radar Topography Mission (SRTM) a DEM with a ground resolution of 90 x 90 m is globally available for free. This resolution is, however, not sufficient for our application. This paper focuses on the generation of a DEM with IBIS-L which is the basis for the combined displacement monitoring of IBIS-L and GPS.

IBIS-L

IBIS-L (Image

         (1)
where topo  is the phase difference due to topography, disp  the displacement in line-of-sight (LOS), atm  the atmospheric disturbance and noise  the instrumental noise and unmodelled effects (i.e. variations of dielectrical properties of the surface).

In contrast to spaceborne SAR, the spatial baseline (i.e. the distance between the two antenna positions in one interferogram) can be selected consciously, depending on the application.

The zero-baseline condition (i.e. spatial baseline is zero) is the ideal condition for displacement monitoring because then topo  is zero.

If the spatial baseline is not zero, a DEM can be derived from topo  . Assuming a spatial baseline only in vertical direction, the relation between elevation z and phase difference topo  of one particular pixel can be approximated by [START_REF] Noferini | DEM by Ground-Based SAR Interferometry[END_REF]

) B r z topo   4  (2)
with  being the wavelength, r the distance to the target and B the vertical baseline.

In case of IBIS-L a spatial baseline length between 0.025 and 0.1 m is possible.

IBIS-L delivers amplitude and phase organized in a two dimensional array of pixels with range and cross-range resolution. Assuming that the position and the orientation of IBIS-L is known, a DEM of the illuminated area is necessary to transform the array of amplitude and phase into the common reference system. The error of the projected positioning x  in LOS with reference to IBIS-L (see Figure 1)

   tan z x  (3)

DEM generation

To determine a DEM with IBIS-L, the phase difference topo  has to be extracted from the observed interferometric phase obs  . As stated in formula (1), additionally to the topography, the observed interferometric phase is also subjected to LOS displacements, instrumental noise and atmospheric disturbances.

Due to the short temporal baseline (i.e. time interval between the images) of about 5 to 10 minutes, the LOS displacements are assumed to be zero. Assuming the instrumental noise to be white noise, the noise can be reduced considerably with the stacking of repeated measurements.

Even with short temporal baselines, the atmospheric effect cannot be neglected. Different methods exist to account for the atmospheric path delay (e.g. [START_REF] Zebker | Atmospheric effects in interferometric synthetic aperture radar surface deformation and topographic maps[END_REF][START_REF] Luzi | Ground-based radar interferometry for landslides monitoring: Atmospheric and instrumental decorrelation sources on experimental data[END_REF].

Assuming constant atmospheric conditions along the path, the path delay r  can, according to [START_REF] Zebker | Atmospheric effects in interferometric synthetic aperture radar surface deformation and topographic maps[END_REF], be approximated by

r T e T P r                 2 1 5 10 73 . 3 10 76 . 7 ( 4
)
whereas P is the atmospheric pressure in hPascal, T the temperature in Kelvin, e the partial pressure of water vapour in hPascal and r the path length.

The assumption of constant atmospheric conditions along the path results in a linear model

r t c r    ) ( (5)
with ) (t c being constant for all pixels at time t and determined as stated in equation (4).

A c c e p t e d M a n u s c r i p t 7 Whether the approach with a linear model is sufficient or not depends on the conditions of the site. Especially the maximum range and the look angle play an important role in the decision of the method for the estimation of atmospheric effects. With great distances and steep look angles non-linearities (e.g. spatial variations of the atmosphere) have to be taken into account.

The corrected phase difference is still wrapped within the range -π to +π. To retrieve the unwrapped phase difference uw  of one particular pixel, the absolute offset (i.e. an integer multiple of 2π) has to be determined

      n uw 2 (6)
with n being the unknown integer.

Here, a conventional branch cut algorithm as described in [START_REF] Goldstein | Satellite radar interferometry: twodimensional phase unwrapping[END_REF] was used for the phase unwrapping. The branch cut algorithm is a path-following method, which works very well for data with high coherence. With a larger spatial baseline a better accuracy can be achieved but the phase unwrapping is more difficult since the height of ambiguity is smaller.

The height of ambiguity is the height z which corresponds to one fringe, i.e. a phase change of   2  topo in equation ( 2) (e.g. [START_REF] Bamler | Synthetic aperture radar interferometry[END_REF]Hartl 1998, Noferini et al. 2008).

DEM generation at a test site

The selected test site is an active quarry in Dieburg owned by the Odenwälder Hartstein Industrie (see Figure 2 a). The coherence of this scene is shown in Figure 2 b).

For the validation of the DEM derived from IBIS-L, a reference DEM was determined by a combination of terrestrial photogrammetry and laser scanning. The laser scan was performed with a Leica HDS 3000 and the photos were taken with the Wild camera P31. The specifications of the used camera and laser scanner are listed in Table 1 andTable 2, 3.

Before the generation of the DEM, the atmospheric effects were accounted for by applying equation (4). A weather station logged humidity, temperature and pressure at the IBIS-L station. Estimations showed that a linear model is sufficient in this test area with a maximum distance of 300 m and an almost horizontal look angle. At this distance the error caused by a change in humidity of 1 % at 25°C is almost 2.5 m in height. Here, during the measurement campaign with duration of 4.5 hours, the temperature rose by about 10°C while humidity dropped by about 20%. This shows the importance of considering atmospheric effects.

Figure 3 shows the observed phase of two control points before and after atmospheric correction. Control point 1 is a triangular reflector (edge length 0.22 m) at a distance of 100 m, control point 2 is quadrangular reflector (edge length 0.22 m) at a distance of 300 m (see also Figure 2). The jumps in the time series are induced by the introduction of a spatial baseline. Since control point 1 is 1 m lower than the IBIS instrument the jump is negative.

Especially in the last section where the observed phase should be zero (because spatial baseline = 0 m), the amount of atmospheric disturbance is clearly visible. The lower graph shows the applied atmospheric correction, whereas the correction for control point 2 is three times higher than for control point 1 (due to the linear model).

The standard deviation of the corrected observation is 0.04 rad for control point 1 and 0.12 rad for control point 2. This leads to an error in height of 0.05 m for control point 1 and 0.11 m for control point 2 with a baseline of 0.1 m. The error doubles for a baseline of 0.05 m.

All possible interferograms were formed and averaged, which resulted in one mean interferogram for each baseline (0.05 and 0.1 m). The DEMs were computed using equation To better understand the differences, the test site is classified into four zones A-D (see Figure 2). Table 4 shows a summary of the different zones and their accuracy. Zone A is a slope with an angle of 40°, which mainly consists of rocks, debris and low vegetation (e.g. grass). The mean coherence is 0.8 ± 0.06 and the found standard deviation 0.8 m. Applying equation ( 3), the positioning accuracy of georeferencing for an elevation angle of 20° is 0.3 m, which is below the spatial resolution of IBIS-L.

Zone B is an old mining area consisting of four terraces with almost vertical walls. The top of each terrace is covered with low vegetation, e.g. bushes and low trees. The mean coherence of this area is with 0.7 not bad but the variation (coherence standard deviation 0.23) is very high due to the covering with vegetation. The standard deviation of height amounts to 2.0 m which leads to positioning accuracy of 0.7 m for an elevation angle of 20°. In general, vertical walls are not represented well in a conventional DEM.

Zone C is mainly covered with vegetation as bushes and trees. The mean coherence is 0.5 ± 0.22. Large discrepancies can be found here, which is expected because vegetation is moving constantly due to wind and thus increases the noise. In this zone the standard deviation amounts to 3.0 m which corresponds to a positioning accuracy of 1.0 m for an elevation angle of 20°.

Zone D only consists of rocks and the coherence is very high but this area was dug off between the measurements of the IBIS-L DEM and reference DEM because more than one A c c e p t e d M a n u s c r i p t 10 month passed by. So, differences of more than ten meters can be found between both DEMs, which are no indication for the measurement accuracy.

Conclusion and outlook

A DEM was generated out of the IBIS-L measurements after the atmospheric effects were removed. A comparison with a reference DEM derived from terrestrial laserscanning and photogrammetry resulted in standard deviations of 0.8 m to 3.0 m depending on the character of the target. The best agreement between both DEMs was found at a slope of rocks and debris. As expected, low coherence targets (e.g. vegetation) decreases the accuracy of the IBIS-DEM considerably. The accuracy of photogrammetry and laserscanning may, however, suffer similar problems in areas with vegetation.

The advantage of using a DEM generated by IBIS-L for georeferencing versus a DEM retrieved from other techniques is clearly that no further measurement instrument is necessary and that the DEM has the same pixel geometry and resolution as the displacement images (i.e. interpolation is unnecessary). The disadvantage is the lower accuracy, especially in areas with low coherence. In areas without vegetation, the positioning accuracy of georeferencing is between 0.3 and 0.7 m in a maximum range of 300 m. Since the maximum achievable resolution in range is 0.75 m with IBIS-L, the accuracy is sufficient.

Further investigations have to be made concerning the accuracy in greater distances. [START_REF] Noferini | DEM by Ground-Based SAR Interferometry[END_REF] achieved a standard deviation of 5.5 m with C-band SAR in 2 km distance, which corresponds to a positioning accuracy of 2 m. As Ku-band SAR, the measurement accuracy of IBIS-L is higher but it also suffers more from atmospheric disturbance and decorrelation.

Here we have shown for the first time the generation of a DEM with a Ku-band GB-SAR IBIS-L and compared it with a high accuracy reference DEM by terrestrial photogrammetry and laser scanning. The goal of the volcano monitoring system is to combine the areal LOS displacements of IBIS-L with 3D displacements measured by GPS. The possibility to Tables Table 1:

Figures

  deformation in real-time with a spatial resolution of several meters and a temporal resolution of several minutes.

  By Interferometric Survey) is a GB-SAR device developed by IDS, Pisa, Italy. It consists of a Stepped-Frequency Continuous-Wave radar unit, operating at 17.2 GHz (Kuband) with a bandwidth of 200 MHz. The synthetic aperture is realized by the movement of the radar unit along a linear rail of 2 m length. The result of one scan along the entire rail is a two dimensional image containing amplitude and phase with a maximum resolution of 0.75 m A c c e p t e d M a n u s c r i p t 5 in range and 4.5 mrad in cross-range (i.e. 4.5 m at 1 km distance to the sensor). The maximum observation range is 4 km. An interferogram is formed by computing the difference between the phase images of two scans at different times. The observed interferometric phase obs  of each pixel is a

  depends on the error of elevation z  and the elevation angle  of the target

  respectively. Both methods delivered point clouds which were registered and combined by transformation with a set of identical control points. The comparison of both DEMs yields a standard A c c e p t e d M a n u s c r i p t 8 deviation of 0.07 m in areas where rock and debris dominate and a standard deviation of 0.41 m where vegetation is present. The combined DEM is shown in Figure 4 a. With IBIS-L, 45 interferograms where acquired with vertical baselines of 0.05 and 0.1 m. The specifications are given in Table

  Figure 4 b. The difference between the reference DEM and the IBIS-DEM is shown Figure 5 a. The mean absolute deviation is 1.8 m. 50 % of the differences are below 0.4 m. Figure 5 b shows the histogram of differences. The corresponding standard deviation is 2.5 m.

Figure 1 :

 1 Figure 1: Influence of the error of elevation on the projected position.

Figure 2 :

 2 Figure 2: Picture (a) and coherence (b) of the observed scene covered by the antenna beam in a quarry in Dieburg. The continuous line is the area of the DEM derived from IBIS; the dashed line represents the area of the DEM derived from photogrammetry and laserscanning.

  2 indicate two control points. The different zones, named A to D are described in the text.

Figure 3 :

 3 Figure3: Phase of the two control points 1 (+) and 2 (x) in 100 and 300 m distance, respectively. The upper graph shows the observed phase measurements (grey) and the phase after atmospheric correction (black). In the lower graph the applied atmospheric correction for both control points is plotted.

Figure 4 :

 4 Figure 4: DEM generated by photogrammetry and laserscanning (a) and by IBIS-L (b).

Figure 5 :

 5 Figure 5: Comparison of DEM derived from IBIS-L and the combination of photogrammetry and laserscanning in m (a). The histogram shows the distribution of differences (b).

  Figure 1

  Figure 2b

  Figure 3

Table 2 :

 2 Specifications of the camera at the test site in Dieburg. Specifications of the laserscanner at the test site in Dieburg.

	Camera	Wild P31
	Focal length	20 cm
	Scanned pixel size	15 μm
	Number of reference points	4
	Baseline	30 m
	Maximum distance	200 m
	Accuracy of height measurement	3 cm
	Laserscanner	Leica HDS 3000
	Instrument type	Pulsed
	Maximum distance	300 m
	Maximum field-of-view	360° x 270°A
	ccuracy of single measurement
	-Position	6 mm at 50 m
	-Distance	4 mm at 50 m
	Used resolution	3 cm at 50 m
	Used max. distance	110 m

Table 3 :

 3 Specifications of IBIS-L at the test site in Dieburg.

	Maximum distance	300 m
	Scan length	2 m

Table 4 :

 4 Found accuracy of the different zones.

	Zone Description

Acknowledgments

The authors thank the Odenwälder Hartstein Industrie for permitting us the installation of IBIS-L in the quarry in Dieburg. This project is funded by the Bundesministerium für Bildung und Forschung within the Geotechnologien program.


James Campbell

Axel Nothnagel

Intra-plate Deformation in West-central Europe

Keywords: Geodetic space techniques, GPS time series, crustal deformation

The central part of Europe north of the alpine orogenic belt is generally seen as a relatively stable area of the western tip of the Eurasian plate. Indeed, up to now, no geodetically significant motions have been detected although an active rift systen running roughly in SSE -NNW direction along the Rhine valley could have some effect on the stability of this region. Presently, the increasing accuracy of geodetic point motions should allow the study of small motions at levels down to nearly 0.1 mm/yr. We start our investigation with a closer look at the 'true' accuracy and significance of GPS derived point velocities of permanent stations. We compare and discuss the different levels of formal errors obtained by the three analysis centers considered in this study (EPN, JPL and SOPAC) and present additional ways of assessing the accuracy using the redundancy offered by different independent analyses and multiple systems operating at one site. On the average, all results indicate that a one sigma level of ±0.3 mm/yr can be seen as a conservative estimate for the horizontal accuracy of point motions in central Europe. On the basis of this assumption we find that at present, the actual velocity field as determined by different analysis groups and centers does not show any significant East-West extensional deformation. We do however see a prominent North-South compressional velocity gradient of about 1mm/yr/1000km (1 nanostrain/yr) which could be associated with the Alpine thrust in conjunction with a south-directed horizontal component of the Fennoscandian Glacial Isostatic Adjustment.

Introduction

In the global scenario of plate tectonics, the continent of Europe at the western end of the Eurasian plate bears its share of crustal deformation: geodetic space techniques and GPS data from permanent stations in particular have been showing prominent site motions on the order of several cm/yr in the Eastern Mediterranean (Turkey and the Aegeis) and on the Mid-Atlantic Ridge (Iceland). The increasing quality and quantity of data in combination with improved analysis systems now also allows to resolve smaller motions down to the millimeter-level and below. A prominent example is given by the observed pattern of horizontal displacements in the range of 0.1 to 2.5 mm/yr due to the Glacial Isostatic Adjustment (GIA) in Fennoscandia [START_REF] Scherneck | BIFROST : Observing the Tree-Dimensional deformation of Fennoscandia[END_REF]). If we narrow our focus to the central area of western Europe, a cursory look at the velocity field does not seem to show any significant motional trends [START_REF] Altamimi | Dense European velocity field and ETRS89 positions and velocities of the EPN stations[END_REF]: central parts appear more or less stable. Geological evidence, however, does present seismo-tectonic activity, in particular on the Cenozoic Rift System which constitutes a zone of crustal weakness extending from the Alps along the Rhine towards the North Sea (Fig. 1). Geologically inferred deformation rates A c c e p t e d M a n u s c r i p t 2 range from average values over about 20 Ma of 0.005 mm/yr to about 1 mm/yr in shorter active phases [START_REF] Ahorner | Erdbeben und jüngste Tektonik im Braunkohlenrevier der Niederreinischen Bucht[END_REF][START_REF] Meghraoui | Active faulting and paleoseismology along the Bree fault, lower Rhine graben, Belgium[END_REF]. Some earlier studies have reported signs of relative point motion in west-central Europe at the level of several tenths of mm/yr, but -as stated by the authors -the data which covered 4 to 6 years of daily GPS observations in the early years of the IGS (1993IGS ( -2000) ) did not seem to be reliable enough to firmly support such claims [START_REF] Nocquet | intraplate deformation in western Europe deduced from an analysis of the International Terrestrial Reference Frame 1997 (ITRF97) velocity field[END_REF]. In their follow-up paper [START_REF] Nocquet | Crustal velocity field of western Europe from permanent GPS array solutions, 1996-2001[END_REF] they state that a level of 0.6 mm/yr (1σ) would approximately draw the limit of significance for site velocities in this area. Similarly, the results of the BIFROST GPS campaign in Scandinavia as published in 2002 can be cited to document a realistic detection limit for GPS derived velocities at that time: the agreement between the observed and modeled motions calculated from tables 2 and 3 in [START_REF] Scherneck | BIFROST : Observing the Tree-Dimensional deformation of Fennoscandia[END_REF] in terms of RMS is σ(v east ) = ± 0.518 and σ(v north ) = ± 0.623 mm/yr. These RMS deviations include the entire GPS error budget as well as GIA model imperfections and local site motions. Considering more recent results on a smaller area of southern Sweden, a set of 12 homogeneously equipped GPS stations with 8 years of 'good' data has produced residual velocity RMS values of σ(v east ) = ± 0.11 and σ(v north ) = ± 0.14 mm/yr [START_REF] Lidberg | An improved and extended GPS-derived 3D velocity field of the glacial isostatic adjustment (GIA) in Fennoscandia[END_REF]). This may certainly be seen as the optimum achievable today.

At present, many of the long standing IGS reference stations as well as some stations from other networks have observational records of ten to fifteen years. The increasing length of time series should have significantly reduced the effects of systematic errors and periodicities in the data. Therefore, it seems appropriate to reconsider the available data and provide an updated picture of crustal motions in the central part of western Europe. In Tab 1 we show a list of central European stations that have been selected for the present study.

We like to add that another motivation for specifically studying the central western part in this publication is due to the fact that we have been carrying out regular campaign style GPS measurements in a section of the Lower Rhine Embayment and the adjacent Rhenish Massif [START_REF] Campbell | Recent movement pattern of the Lower Rhine Embayment from tilt, gravity and GPS data[END_REF]. This project has been primarily aimed at monitoring mining subsidence but also includes the investigation of tectonic motions connected with the Eifel-Ardennes uplift and horizontal extension on the northern part of the Cenozoic Rift System. First results from 8 annual campaigns (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) show a small SW -NE extensional motion for two pairs of stations in the Rhenish Massif across the Lower Rhine Embayment of ~0.5 ± 0.5 mm/yr, a possible trend that is still buried in the noise. The annual campaigns are continuing in an extended net of sites providing results from a time span of more than 15 years [START_REF] Görres | Recent Site Motions in the Lower Rhine Embayment and the Eifel from 15 years of GPS data[END_REF].

Accuracy of GPS derived site velocities

At the site velocity level of 1 mm/yr and below, the well known systematic errors endemic in geodetic Global Navigation Satellite Systems (GNSS) are causing growing problems. Without going into too much detail we list the most important error types [START_REF] Zhang | Southern California permanent GPS geodetic array: error analysis of daily position estimates and site velocities[END_REF][START_REF] Nikolaidis | Observation of geodetic and seismic deformation with the Global Positioning System[END_REF][START_REF] Steigenberger | Reprocessing of a global GPS network[END_REF], noting that the changes in time of these sources are most harmful: reference frame realizations (changes of), satellite orbits (changes, improvements), atmospheric delay (new models), antenna calibration (changes, new standards), multipath environment (changes in antenna setup, neighboring structures). If we focus our attention on the detection of crustal motion we also have to consider ground movements of local origin as well as regional effects such as atmospheric loading and hydrological variations.

A c c e p t e d M a n u s c r i p t 3 While research on each one of these error sources is going on, it is also appropriate to make use of the redundancy and variety of analyses from different groups and centers, the results of which are often publicly available. Certainly the most important outflow of these activities may be seen in the realization of the International Terrestrial Reference Frame (ITRF) which combines all the information contained in the individual results into one consistent global product (Altamimi et al. 2007). The ITRF includes site cordinates at a reference epoch as well as site velocities together with their respective errors (variance/covariance information).

Normally, it would be sufficient to use these errors and discuss the listed site velocities in terms of their geotectonic significance. If we are faced, however, with a rather small subset of sites that seem to be moving randomly at rates between 0.1 and 1 mm/yr, it is necessary to examine more closely each one of these sites and try to establish the quality and consistency of both the data and the subsequent analysis. To a certain degree, this has of course been done in the compilation of the successive ITRF realizations, although some particular aspects or problems of the individual analyses require closer examination in deeper studies. Thus, for the present investigation, we examine several separate solutions before the global combination takes place. For ease of data access we chose results from GPS analyses of three different organisations, i.e. the European Permanent Network (EUREF EPN), the JPL Flinn analysis and the SOPAC data center (Tab. 1a). These results have a high degree of independency because they are produced with entirely different software systems and follow their own particular analysis strategy. In addition, we will briefly address recent results of analyses of European VLBI observations [START_REF] Campbell | European VLBI for crustal dynamics[END_REF] carried out at our institute (IGG). With this material we will try several different approaches to provide additional information about the accuracy of the site velocities. We consider comparisons between results from different independent analyses, the temporal evolution of results at the centers as data accumulate and software systems are improved, and use the results from different receivers at one site. To start with, we need to take a closer look at problems related to the application of the stochastic models commonly used in the GPS analyses.

Formal errors

Generally, station coordinates and velocities (rates of change of coordinates) are computed simultaneously in so-called global solutions together with their respective variances and covariances (e.g. [START_REF] Nikolaidis | Observation of geodetic and seismic deformation with the Global Positioning System[END_REF]). In the case where site velocities are derived directly from the time series of daily (or weekly) coordinate solutions, the velocity errors come as a byproduct from the least squares model fits to the data. Because in this process the data covariance information is not propagated, it has to be reproduced empirically from the time series itself. Neglecting the correlation between data points by setting cov(τ) = 0 would mean to stipulate the domination of white noise. In this case, for a series of N equidistant data points in time the root variance ( ν σ ) of the rate in a two-parameter linear fit becomes [START_REF] Zhang | Southern California permanent GPS geodetic array: error analysis of daily position estimates and site velocities[END_REF]): Facing the systematic effects seen in the majority of the time series, analysis groups have developed more refined stochastic models based on covariance input derived from spectral analysis of the series [START_REF] Zhang | Southern California permanent GPS geodetic array: error analysis of daily position estimates and site velocities[END_REF][START_REF] Nikolaidis | Observation of geodetic and seismic deformation with the Global Positioning System[END_REF][START_REF] Williams | The effect of coloured noise on the uncertainties of rates estimated from geodetic time series[END_REF]. This leads indeed to more realistic error estimates, although the determination of parameters for the different types of noise remains rather uncertain [START_REF] Mao | Noise in GPS coordinate time series[END_REF]).

) t / T 2 )( t / T 1 ( t / T 12 T r ∆ + ∆ + ∆ σ = σ ν (1)

Power law noise

Spectral analysis is a useful tool for analysing the noise features typically seen in most of the time series [START_REF] Agnew | The time-domain behavior of power-law noises[END_REF]). On the basis of a generalised power law of the form P = P 0 (f/f 0 ) κ with f as sampling frequency and P 0 and f 0 scaling coefficients, three types of noise can be defined, characterized by values of 0, -1 and -2 of the exponent κ, which stand for white noise, flicker noise and random walk, respectively. Introducing covariance information estimated from the power specra [START_REF] Zhang | Southern California permanent GPS geodetic array: error analysis of daily position estimates and site velocities[END_REF], the formal errors obtained in the fit increase with growing values of κ , which is easily appreciated knowing that data correlation reduces the redundancy of a least squares solution. To first approximation and for n equidistant data points the rate errors are proportional to 1/n√n for white noise, 1/n for flicker noise and 1/√n for random walk.

A comparison of the velocity errors from three different analysis centers in the histograms in Fig. 3a,b shows the effect of applying different stochastic models: the JPL Flinn Analysis apparently uses a white noise model while EPN and SOPAC use coloured noise models.

As an example of spectral analysis, we show the power spectrum of the JPL time series of station Wettzell (WTZR), displaying the annual peak as well as peaks at higher frequencies such as the fortnightly tidal oscillation of 13.66 d, the latter indicating a deficiency in the tidal model of the analysis software (Fig. 4).

Here we like to point out that the high frequency oscillations are of little or no effect on the rate estimates of longer series, quite in contrast to the detrimental effect of the lower order harmonics (Blewitt and Lavallee 2002). To document this, in a special low order Fourier analysis we chose the base frequency ω 0 = 2π/T 0 with T 0 = 20yr for 10-yr time series (Fig. 5a,b) obviating the large energy in the low period range from 1 to 20 years. As an example, we choose the BRUS North component of the 13,4 yr JPL data set after velocity estimation showing the detrended residuals (Fig. 5a). The presence of long period signals is obvious and could be related to a low frequency type of oscillation. The nature of periodic or quasiperiodic signals beyond the ubiquitous annual oscillation still remains elusive. It could just be the result of changes in the equipment and/or analysis approach. Other sources could be A c c e p t e d M a n u s c r i p t 6 found in anomalous hydrological cycles (e.g. groundwater variations) or monument instability. If groups of stations suffer similar variarions one should consider the impact of higher order ionospheric terms [START_REF] Steigenberger | Reprocessing of a global GPS network[END_REF]. Irrespective of the origin of the longer period variations one can apply the rules discussed by [START_REF] Blewitt | Effect of annual signals on geodetic velocity[END_REF] to estimate the velocity bias resulting from the unmodelled oscillation. Thus, in the case of a 10yr 'wave' with an amplitude of 2.5 mm one would obtain a maximum bias of 0.5 mm/yr if the data span T matches the length of the period, i.e. 10 yrs, but the bias would quickly fall to zero for T = 15 years, oscillating around decreasing values from 0.25 mm/yr downwards for longer T (scaling Fig. 3 in Blewitt and Lavallée).

Problems with individual stations

The assessment of the quality of fit of the time series of the individual analysis centers is made easy by the steadily improved content and layout of the websites (Tab.1). Moreover, it is possible to download the data and apply one's own models and software to reproduce and check the results (Section 3.2).

Going through the time series plots offered by the different websites, we found several aberrant results, most of which were corrected in the ensuing updates. Here, we like to show just one example which is typical for any one of the analysis centers and a consequence of the enormous amount of data that has to be routinely treated in global networks. In the SOPAC analysis of the station Kootwijk (KOSG) the data of the GIG campaign in 1991 and the data in early 1992 caused difficulties in the fitting model that led to velocity estimates off by more than 1 mm/yr. The problem could be easily identified in the de-trended residual plots (Fig. 6) and was subsequently corrected by excluding the data from 1991 to 1992.3. This, however, was the only correction actually applied, in order to safeguard the authenticity of the data drawn from external sources.

Temporal evolution of velocity estimates

To evaluate the accuracy of the velocity solutions, it is instructive to look at the effects of improving data quality, software upgrades, such as new and better models for atmospheric delays, antenna calibration on the results as they evolve. Here we like to present as an example several different downloads from the JPL FLINN Analysis Center that occurred during our work in the period 2000 to 2008 (Fig 7a,b). From our list (Tab. 1) we selected seven stations in central Europe with best performance in terms of length of series and quality of velocity fit.

The salient features in the behaviour of the seven selected stations are common shifts related to changes in the reference frame realised in the different solutions. For studies on crustal deformation in regional and local networks these changes are of little or no effect. The important aspect here is the change in the relative behaviour of the velocity values among the stations. The trend observed in Fig. 7a,b is evident: most velocity values run in parallel and show relative variations smaller than 0.3 mm/yr. To better quantify our conclusion, we formed differences between pairs of stations to eliminate the frame effect and calculated a WRMS from differences between the two latest download dates of 2007. From this example we conclude that one-sigma velocity errors cannot presently be considered to be smaller than 0.1 mm/yr for the best stations. Of course, larger variations are seen in the first sections between 2000 and 2003, in compliance with the fact that the time series were still rather short (4 to 6 years) in the early phase of permanent observations.

Different receivers/antennas at one site

Operating different systems at one site has proven to be a very efficient way of establishing realistic accuracy estimates. Within the IGS network the Geodetic Observatory of Wettzell offers the possibility to compare velocity results from five different receiver/antenna assemblies deployed on a platform of the observing tower of the main building (Söhne et al. 2005, Fig. 8a). The data from all five of these receivers are being routinely processed and made publicly available at the JPL Flinn analysis center (ref. see Tab.1). In spite of the small sample size, we believe that the comparison does provide useful complementary information for the present accuracy discussion.

The errors that can be brought to light by comparisons at one single site are mainly due to the instrumentation and include antenna calibration problems, receiver malfunctions or other deficiencies. If we focus on the velocity determination and not so much on the coordinates themselves, it is the change in any of the system parameters that is essential. One can imagine changes in the antenna setup, the mount and the environment as being most harmful to site velocity determinations.

In Fig. 8b we plotted the velocity vectors for the five receiver/antenna assemblies with reference to the weighted mean of the individual values in order to illustrate the spread of the

Comparison between solutions from different GPS analysis groups and centers

In the comparisons, we have narrowed our focus on the West-central part of Europe and plotted the vectors in a magnified scale in order to be able to identify any obvious streams of velocities (Figs 9a,b,c). The global reference frames of the JPL and SOPAC velocity data were rotated to the Eurasia-fixed frame using the absolute plate rotation values of the ITRF2005 given by Altamimi et al. 2007. In order to facilitate the visual comparison between solutions, we had to apply a small additional rigid rotation to the SOPAC data set. For the purpose of this study, i.e. comparing only deformations or relative site motions, the actual amount of rotation is irrelevant. Moreover, more sophisticated approaches using adjustments on a selected set of 'stable' sites would require a significantly larger set of points (Scherneck et al. 2002 , Nocquet and[START_REF] Nocquet | Crustal velocity field of western Europe from permanent GPS array solutions, 1996-2001[END_REF].

In all three representations (EPN, JPL and SOPAC), the E-W pattern of motions remains rather accidental and does not suggest any significant E-W-extension above the level of about 0.3 mm/yr. In contrast to this we see a rather persistent north-eastward trend of the sites in the northern Alps and their foreland. In southern Sweden, the site of Onsala, whose motion is also well documented in the SWEPOS and BIFROST campaigns, shows a sustained south-western trend. This motion is in good agreement with GIA models [START_REF] Scherneck | BIFROST : Observing the Tree-Dimensional deformation of Fennoscandia[END_REF]. The stations in central Germany (WTZR, POTS, HOBU, PTBB) all have rather small residual N-S In this context we should like to emphasise that we have to refrain from any deeper discussion of the vector fields of motion, because they are highly dependent on the choice of a 'crust fixed' reference frame. [START_REF] Nocquet | Crustal velocity field of western Europe from permanent GPS array solutions, 1996-2001[END_REF] went through great pains to define a stable 'no motion' reference in central Europe using all available stations and had to accept a remaining level of residual motion of 0.5 mm/yr. It is of interest to note that even at this early stage (data from 1996-2001) the coherent N-E motion of the three pre-Alpine sites ZIMM, HFLK, GRAZ and the N-S convergence is just about to manifest itself (their Fig. 3).

To arrive at quantitative measures that can be used in further comparisons, we computed mean velocity gradients from linear fits in E-W and N-S profiles, or rather profile swaths for the available stations in each one of the analyses (Figs 10a,b,c), in a way similar to e.g. [START_REF] Mourabit | GPS-Geodynamic research in and adjacent to Morocco[END_REF]. The deformation trends are quite consistent in all three of the solutions (Tab. 4) and support the visual impression from the vector plots in Figs. 10-13, i.e. that the E-W extension is all but marginally recognisable, while the N-S compression appears quite clearly with a mean value over 1000 km of 1.2 ± 0.15 mm/yr.

Tab. 4: Average strain rates in N -S and E -W profiles. The gradient units are equivalent to strain rate in nanostrain/yr, i.e. 10 -9 /yr. n = number of stations in each data set. Mean = weighted mean of all three data sets with n used as weight.

Comparison with VLBI measurements

As far as independent techniques such as SLR and VLBI are concerned, the small size of the area under study undoubtedly precludes any useful comparisons with the GPS results. In the frame of the regular European VLBI sessions, however, there is a possibility to look at a three station configuration which covers the central part of western Europe, i.e. the triangle formed by the stations of Onsala, Wettzell and Effelsberg (Fig. 11). The Effelsberg 100m radio telescope is located just west of the Rhenish Embayment in the northern foothills of the Eifel-Ardennes Massif. Although heavily burdened by astronomical observations, regular time windows have been set aside for geodesy to permit a meaningful contribution to the global and European geodetic VLBI programs. The long time series in combination with the low level of long period systematic effects makes the VLBI results particularly suitable and valuable for the present study.

In particular, the time series for the length of the baseline Effelsberg -Wettzell (Fig. 12), which crosses the Rhenish rift zone, gives an indication of the VBI data rate and quality and leaves little room for any E-W-extensional motion beyond the level of 0.3 mm/yr (in terms of integrated relative motion over 500 km). From the vectors in Fig. 10 we can derive the

Conclusions

The overall picture gleaned from the comparisons may be interpreted as a consequence of the presence of a N-S-compressional regime due to the alpine thrust in the South and the horizontal component of the GIA in the North. We found signs of an N-S velocity gradient of 1.2 mm / 1000km = 1.2 nanostrain per year, or in other words a compressional deformation at a level of 1 to 2 mm/yr between southern Scandinavia and the northern Alps and their forelands. This observed trend is inherent in both the GPS analyses as well as in the VLBI results. Furthermore, we can confirm earlier research on the absence of any detectable E-W extensional motion [START_REF] Nocquet | Crustal velocity field of western Europe from permanent GPS array solutions, 1996-2001[END_REF], and lower the limit of significance to a value of ≤0.3 mm/yr. This is an important result in view of the seismo-tectonic risk assessment in the area of the Cenozoic Rift System [START_REF] Camelbeek | Large earthquakes in Northern Europe more likely than once thought[END_REF]. However, to guard against any misconceptions about the interpretation of the observed discrete site motions we emphasise that the results stated here represent integral relative motions over large distances and should not be compared directly with locally derived motions on fault systems.

The group of stations in the north-west of central Europe shows several larger particular motions that do not seem to fit with any general tectonic trends. This behaviour forms a distinct contrast to the stations further east, which have been showing a rather stable relationship for many years now.

The small number of sites involved in this study still asks for caution before any firm conclusions can be drawn, in particular considering the unknown influence of long period effects discussed in section 3.2. Therefore, the inclusion of more permanent stations in the crustal deformation studies is urgently needed. In this respect, the German GREF network will be of great benefit [START_REF] Ihde | National Report of Germany. Mitt. des Bundesamtes für Kartographie und Geodäsie[END_REF], much like the French Regal network and others to be used in combinations for creating denser fields of points for high resolution velocity determination [START_REF] Dong | Estimating regional deformation from a combination of space and terrestrial geodetic data[END_REF][START_REF] Nocquet | Crustal velocity field of western Europe from permanent GPS array solutions, 1996-2001[END_REF][START_REF] Camelbeek | Crustal deformation in stable continental Europe: a comparison of seismicity, geodetic and geologic information[END_REF].

  of the series in years, ∆t data interval and r σ the error of a daily coordinate solution. In order to account for the effect of grossly underestimating the actual error level, the ν σ -values may just be scaled by an empirical factor which can be determined e.g. from test campaigns. In Fig.2have we used a scale factor of 3.0 to be applied to the typical white noise error value of ±3.0 mm of a single day GPS horizontal position determination. In this representation of the basic error law for time series of length T(y) = N/365 we have included the actual errors from the EPN (CLEAN) time series analysis downloaded in 4/2008 for a subset of 38 stations. Those stations closest to the line of the upscaled error law represent the 'best behaved', while those further away from the ideal error line are more strongly affected by anomalous signals, a situation that can be readily verified by examining the corresponding time series plots.

  4 and 2008.3 (Tab. 2a,b).

  Evolution of velocity results of the JPL analyses. DD = Double differences for the two latest download dates T 1 and T 2 . The resulting WRMS for the North and East components are: wrms North = ± 0.101 and wrms East = ± 0.113 mm/yr. σ T1 and σ T2 are the white noise rate errors from the respective

  the formal errors (see section 3.1) were correct, statistically two of three vectors (65%) should vary within a circle of only 0.1 mm radius. The actual spread reflects above all the effect of observing time, i.e. different length of the series. Stations with series of less than 5 years suffer from anomalous deviations of more than 1 mm/yr(Blewitt and Lavallee 2002). In Tab. 3a,b we compiled the RMS errors for the north, east and up components, showing that a standard time series of 10 years duration would realistically produce RMS errors of about 0.3 mm/yr in the horizontal components.

  sites west of the Rhenish rift system in Belgium, the Netherlands and southern Britain (HERS) tend to show somewhat larger individual motions although they are located fairly closely together.

  the E-W-direction, i.e. dv E-W ≤ 0.3 mm/yr. Conversely, we do see part of the N-S compression observed by GPS, i.e. about 1 mm/yr.

Fig. 9a :Fig 11 :

 9a11 Fig. 9a: EPN velocity vectors for central-western European stations; shaded lines indicate the limits of the E-W and N-S profiles

FigureFig. 1 :

 1 Figure

Figure

  Figure

Fig 2 :Fig 4 :Fig. 5a :Fig. 5b :

 245a5b Fig 2: Time series error law for white noise (enhanced by an empirical factor 3) for comparison with errors given in the EPN CLEAN velocity results.

FigureFig. 6 :

 6 Figure

FigureFig 7a :

 7a Figure

Figure

  Figure

Figure

  Figure

Figure

  Figure

Fig. 10a :

 10a Fig. 10a: N-S profile and E-W profile velocity gradients from EPN data

Figure

  Figure

Fig

  Fig 10b: JPL N-S and E-W mean velocity gradients

Figure

  Figure

Fig. 10c :

 10c Fig. 10c: SOPAC N-S and E-W mean velocity gradients

Figure

  Figure

Fig. 12 :

 12 Fig. 12: VLBI baseline lengths, Effelsberg -Wettzell, slope = -0.086 ± 0.12 mm/y

Figure

  Figure

Tab. 1a: Data sources.

	M a n u s c r i p t
	Data Source EPN -GPS JPL -GPS SOPAC-GPS IGG-VLBI Tab. 1b: Stations used in this study Institution EUREF Permanent Network www.epncb.oma.be Access JPL Flinn Analysis sideshow.jpl.nasa.gov Scripps Orbit and Permanent Array Center sopac.ucsd.edu Institut für Geodäsie und Geoinformation, Bonn A c c e p t e d	Analysis software Bernese GPS GIPSY GAMIT Calc/Solve
	4

Acknowledgements.

We appreciate the opportunities to use the 100m telescope of the Max Planck Institute for Radio Astronomy at Effelsberg in the European geodetic observing sessions. The VLBI data analyses have been carried out by Thomas Artz (IGG). Moreover, we like to acknowledge the availability of data and velocity solutions by the analysis centers of the European Permanent Network, the JPL Flinn Analysis and the Scripps Orbit and Permanent Array Center (SOPAC).