GEOLOGICA CARPATHICA
, DECEMBER 2016, 67, 6, 509 – 524
doi: 10.1515/geoca-2016-0032
www.geologicacarpathica.com
Surface strain rate colour map of the Tatra Mountains
region (Slovakia) based on GNSS data
MARTIN BEDNÁRIK
1
, JURAJ PAPČO
2
, VLADIMÍR POHÁNKA
1
, VLADIMÍR BEZÁK
1
,
IGOR KOHÚT
1
and LADISLAV BRIMICH
1
1
Earth Science Institute, Slovak Academy of Sciences, Dúbravská cesta 9, P.O. Box 106, 840 05 Bratislava, Slovakia; geofmabe@savba.sk
2
Department of Theoretical Geodesy, Faculty of Civil Engineering, Slovak University of Technology, Radlinského 11,
813 68 Bratislava, Slovakia
(Manuscript received September 29, 2015; accepted in revised form September 22, 2016)
Abstract: The surface deformation of the Tatra Mountains region in Western Carpathians can nowadays be studied
directly thanks to precise geodetic measurements using the GNSS. The strain or stress tensor field is, however, a rather
complex “data structure” difficult to present legibly and with sufficient resolution in the form of a classical map. A novel
and promising approach to the solution of this problem is coding the three principal strain or stress values into the three
colour channels (red, green, blue) of an RGB colour. In our previous study, the colour depended on the stress tensor shape
descriptors. In the current study, the adapted colouring scheme uses a subset of shape descriptors common to stress and
strain, which differ only in the scaling factor. In this manner, we generate the colour map of the surface strain rate field,
where the colour of each grid point carries the information about the shape of the strain rate tensor at that point.
The resulting strain rate colour map can be displayed simultaneously with the map of the faults or elevations and be
easily checked for the data or interpolation method errors and incompatibility with the geophysical and geological
expectations.
Keywords: GNSS, surface deformation, tensor field visualization, three-dimensional colour function, Tatra Mountains.
Introduction
Tatra Mountains
The Tatra Mts. (Fig. 1) are the northernmost core mountain
range in the arc of the Western Carpathians. Tectonically, they
belong to the Inner Western Carpathians (recently reaffirmed
by Bezák et al. 2004). Their present-day morphology is
a product of the most recent tectonic processes in the Miocene
to Quaternary. After the oblique collision of the block of the
Inner Carpathians with the European Platform and the subduc-
tion of Magura Ocean reached standstill, new tectonic pro-
cesses ini tiated in the mantle have led to increased volcanic
activity and formation of horst-graben structures, with most
horsts uplifted asymmetrically. The asymmetrical uplift of the
block of the Tatra Mts. proceeded along the Sub-Tatra fault
band (e.g., Nemčok et al. 1994). The Tatra Mts. horst com-
prises in itself the tectonic units from the previous orogenic
periods, mainly the Hercynian and the Palaeoalpine. In the
Hercynian period, crystalline units (metamorphics and grani-
toids) were formed. During the Palaeoalpine Cretaceous
period, the superficial nappes of Mesozoic complexes were
formed, preserved now mainly on the northern side of the
Západné and Vysoké Tatry Mts. and in the Belianske Tatry
Mts. A detailed view of the geology of the Tatra Mts. is given
by Nemčok et al. (1993). The Neoalpine to recent tectonics of
the Tatra Mts. have been studied by, for example, Sperner &
Ratschbacher (1995), Burchart (1972), and recently by
Králiková et al. (2014). A valuable summary of the research
related to the exhumation of the Tatra Mts. can be found in the
introductory chapter of Anczkiewicz et al. (2015).
Deformation measurements in the Western Carpathians and
Pannonian Basin
The completely blank space of Slovakia on the World Stress
Map (WSM) compiled by Heidbach et al. (2009) represents
the reluctance of the local researchers to contribute to the WSM
rather than a total lack of any research. The early history of
stress/strain field research in Western Carpathians and Panno-
nian Basin definitely includes the extensometric measure-
ments. A summary of their results can be found in Mentes
(2008). Their spatial resolution is very low (only five extenso-
metric stations for the whole area of interest) and they only
give the projection of the strain tensor at a site usually into
one, exceptionally two approximately perpendicular horizon-
tal directions. As the extensometric measurements are sampled
frequently enough in time, they could be able to detect sudden
tectonic stress changes at the site. Nevertheless, due to the
environmental influences (change of temperature, preci-
pitation etc.), they never did so with sufficiently high
reliability.
Another category of strain measurement devices frequently
used in Central Europe is the dilatometers (crack gauges),
especially the mechano-optical ones of the type TM-71
(Košťák 1969). These register the relative motion (in all its
510
BEDNÁRIK, PAPČO, POHÁNKA, BEZÁK, KOHÚT and BRIMICH
GEOLOGICA CARPATHICA
, 2016, 67, 6, 509 – 524
degrees of freedom) of two blocks of rock across the crack that
separates them. The favourite installation sites are the under-
ground cavities of caves or mine galleries (e.g., Briestenský
2008; Briestenský & Stemberk 2008), thus giving a rise to
a new branch of active tectonics research (Littva et al. 2015).
A study by Fojtíková et al. (2010) relates the focal mecha-
nisms of the local earthquakes to the stress field in the Malé
Karpaty Mts. Given the number and spatial configuration of
seismometers needed to determine the focal mechanisms of
many local events, and the uneven spatial and temporal distri-
bution of seismic activity in the Western Carpathians, we cannot
hope that any larger portions of Western Carpathians will ever
be surveyed in this manner.
With the onset of new technologies, much effort is exerted
to measure the surface movements directly. The comprehen-
sive dissertation of Papčo (2010) directly applies to the same
geographic region of interest as this paper. The data output of
the study serves as our source data. Our method provides
a further enhancement of their geophysical readability and
interpretability. As Papčo’s dissertation is not widely availa-
ble, in the next chapter, we give a short summary focusing on
information relevant to this paper.
GNSS measurements in the Tatra Mountains region
Global navigation satellite systems (GNSS) can be used as
very efficient tools for monitoring the velocity of lithospheric
plates, including the inner plate motions.
The Tatra Mountains region has
attracted the specialists in geodetic
monitoring since the early 1980-s.
The pioneering work has been
done by Hradilek with his high
precision
classical
geodetic
measu rements (Hradilek 1984).
The modern period of this long-
term activity started in 1997 within
the project CERGOP (Czarnecki
& Mojzeš 2000) and CERGOP2 -
Environment (Fejes 2002) by
establishing the special GNSS
geodynamics network TATRY.
The whole network consists of 30
points, which can be divided into
two main groups: primary geody-
namics points (18 sites), and the
rest of the network dedicated to
other research purposes (e.g., local
gravity field modelling). Geody-
namics points are located in the
highest alpine parts (points ZISE
(Žiarske sedlo), KAWI (Kasprowy
Wierch), PRSO (Predné Solisko)
and LOMS (Lomnický štít)), val-
ley parts (points ZICH (Žiarska
chata), HAGA (Hala Gąsie ni-
cowa), TIDO (Tichá dolina), BIEL (Bielovodská dolina),
SDOM (Sliezsky dom), HREB (Hrebienok), SKPL (Skalnaté
pleso) and TAKO (Tatranská kotlina)), and 6 points are situated
outside the main area (GANO (Gánovce), ROHA (Roháčka,
near Liptovský Mikuláš), LIES (Liesek), KACI (Kacwin),
RWOO (Rolów Wierch) and GUBA (Gubałówka)). Precise
epochal GNSS measurements have been performed annually
since 1998 with the durations of approximately 96 hours.
Processing was done in the Bernese software version 5.0
(Dach et al. 2007) in accordance with the EPN and CEGRN
recommendations (Stangl 2007) using precise IGS orbits, EOP
parameters (IGS 2016), absolute GNSS antenna phase centre
information and Niell mapping function for the tropospheric
delay modelling. All campaigns (from 1998 to 2008) were
processed in the ITRS global coordinate system. Connection to
ITRF2005 (Altamimi et al. 2007) was established by using the
surrounding permanent GNSS stations of EPN: BOR1, GOPE,
GRAZ, JOZE and PENC (EPN 2016).
Final coordinates and velocities of the points (in ITRS
global coordinate system) were estimated in the special pro-
gram COMB_NET, which was designed for combinations of
networks taking into account the time epoch (Hefty 2004;
Hefty & Gerhátová 2006). The inner plate velocities were
derived from the Eurasia lithospheric plate velocity model
APKIM2005D (Drewes 2009).
The horizontal velocities reach up to 2 mm/year, but at the
majority of the points they are below 1 mm/year (cf. Appendix;
Fig. 2). Vertical velocities of the points are very unstable and
Fig. 1. Tectonic scheme of the Tatra Mts. region based on Bezák et al. (2008). The coordinates [m] of
its corners are given in UTM Zone 34 North reference system. In the upper left corner, the overview
map of Slovakia shows the studied area with a black rectangle.
511
SURFACE STRAIN RATE COLOUR MAP OF THE TATRA MOUNTAINS
GEOLOGICA CARPATHICA
, 2016, 67, 6, 509 – 524
.
(1)
The dot over the symbols denotes the first time-derivative of
strain and displacement. The indices after the comma denote
partial derivatives with respect to the l-th or j-th spatial
coordinate.
Of that velocity field
, we have available only a few sam-
ples — the horizontal surface velocities measured at the sta-
tions. From those samples, taking into account their errors, we
have to interpolate the horizontal velocity field covering the
whole region. The interpolated velocity field can be specified
in the form of a function formula, which can be differentiated
analytically, or as a regular grid of values, which can be diffe-
rentiated numerically.
At the stations’s positions, the exact interpolation methods
yield exactly the same values as the respective measured values.
This is not very desirable in the situation of known, relatively
large data errors. The approximate interpolation methods,
which do not adhere to the measured values exactly, should
then be preferred. This was also our choice.
The next choice depends on our notion of regional dyna-
mics. A region with a stable crust can be regarded as moving
as a whole, with small, regional scale perturbations. In this
case, we shall search for a whole region interpolation formula.
In a less quiet region, the local dynamics concentrated to the
vicinity of faults and other active structures may locally prevail
over the regional trends. In the latter case, the local interpola-
tion methods, effectively taking into account only the nearest
stations, may be a better choice. An increase of the (polyno-
mial) order of the whole region interpolation formula could
have a seemingly similar effect of
capturing the local behaviour.
Nevertheless, this solution has
many side effects. The worst of
them is a better adherence to the
values measured at the sites at the
expense of very unsatisfactory
behaviour elsewhere in the region.
Increasing the order of approxi-
mation makes more sense in the
local interpolation methods.
We tested both the whole region
least squares interpolations, im -
ple mented in a script written in
Mathematica 9 by Wolfram Re-
search, and the local moving least
squares interpolation (Pohánka
2005), implemented in a stand-
alone program.
In the whole region weighted
least squares interpolation under
Mathematica 9 we took advan-
tage of the built-in Minimize[ ]
function to find the values of
the coefficients c
mn
of the
poly nomial
(
)
2
,
,
j
l
l
j
l
j
u
u
+
=
e
j
u
problematic, because none of them surpasses its uncertainty
and they show very variable trends. Next analysis and strain
field determination were performed without the points with
problematic behaviour (BIEL, ZICH, TIDO, PRSO, GUBA,
RWOO and HREB). The main reasons for these exclusions
were locations of the points (sites affected by landslide
or moraine movements), character and magnitude of the
velocity, sudden jumps in the time series or shorter measure-
ment period.
For the purposes of this paper, two more stations — GANO
and LIES — were excluded from the dataset based on consi-
derations about the sites’ geological settings. GANO can be
influenced by slope motions. LIES is located on Neogene
sedi ments. TAKO, situated on glaciofluvial deposits, has been
kept in the dataset. The difference is that the contribution of
TAKO to the interpolation function can be corrected by contri-
butions of the neighbouring stations (and it has also a lower
weight in the weighted least squares due to its higher standard
deviation of velocity), whereas LIES is a stand-alone site that
would influence a broader area.
Methods
Interpolation methods
The main task is to produce a colour map of the surface
strain rate tensor field. The strain rate tensor contains the spa-
tial derivatives of the surface velocity field. In Cartesian coor-
dinates and index notation
Fig. 2. Measured surface velocities (arrows) and their standard deviations (ellipses) within the area
of interest with an extent of 70 × 47.5 km
2
(the coordinates [m] of its corners are given in UTM Zone
34 North reference system). Here, the different colours are used only as a graphical tie between the
measurement site and the standard deviation of velocity at the site. In the upper left corner, the over-
view map of Slovakia shows the studied area by black rectangle.
512
BEDNÁRIK, PAPČO, POHÁNKA, BEZÁK, KOHÚT and BRIMICH
GEOLOGICA CARPATHICA
, 2016, 67, 6, 509 – 524
(2)
that yield the minimum of
∑
-
k
O
k
I
k
k
v
v
w
2
, where
(
)
k
k
I
I
k
y
x
v
v
,
=
is the interpolated value and
O
k
v
is the value
observed at the site (x
k
, y
k
) and w
k
, w
k
>0 , are the weights at
the sites. v
I
(x, y ) stands for any of the horizontal surface velo-
city components. In our implementation, the weights of the
approximation errors are proportional to the reciprocals of the
specified standard deviations of the measured velocity values
at the sites. Thus, the smaller the standard deviation, namely
the more reliable the value, the greater weight is attributed to
the site. The site’s weight could also involve a station’s geody-
namical reliability factor based on an expert’s assessment of
the station’s geological setting. For instance, the stations on
unconsolidated sediments could get a lower rating. We did not,
however, implement the station’s rating. We instead excluded
the problematic stations GANO and LIES from the dataset.
The weights based on standard deviations did not have a very
significant influence on the results, as the standard deviations
vary from site to site only moderately (with an exception of
TAKO — cf. Fig. 2) and the almost equal weights did not
favour any of the stations.
Our iterative least squares algorithm with adaptive weights
yielded interesting results. Its basic idea is simple: if the
approximation error
( )
O
k
I
k
v
i
v
-
in the i-th step at the k-th site
exceeds the standard deviation
d
k
(here we deviate from the
usual notation of standard deviation by
s
k
, as this symbol
denotes principal stress value elsewhere in the paper), then in
the i +1-th step the site’s weight w
k
(i + l) should be increased
by the factor of
( )
( )
k
O
k
I
k
k
v
i
v
i
a
d
-
=
(3)
to
( )
( ) ( )
i
w
i
a
i
w
k
k
k
=
+
1
(4)
in order to reach a more favourable a
k
(i + l) ratio (i.e. a lesser
exceedance) in the i + l-th step of weighted least squares
approximation. In this manner, the weights of all sites shall be
iteratively adapted. The iteration process stops when for all
sites
( )
( )
O
k
I
k
O
k
I
k
v
i
v
v
i
v
-
=
-
+
1
.
The observed convergence
of the process was rather slow (up to hundred iterations
needed) and very moderately uneven. It seems to be, however,
independent of the initial weights w
k
(0) . We tried w
k
(0) =1 as
well as w
k
(0) =
d
k
–1
.
Let us denote
( )
i
v
v
I
k
i
A
I
k
∞
→
=
lim
. The adaptive weights algo-
rithm tends to minimize the maximum approximation error
O
k
A
I
k
k
v
v -
max
reached within the dataset, whereas the origi-
nal least squares approximation directly minimizes
∑
-
k
O
k
I
k
k
v
v
w
2
. Interestingly, the Mathematica’s Minimize[ ]
function was unable to find the global minimum of
O
k
I
k
k
v
v -
max
directly. The detour through iterative least
squares is a necessity.
The moving least squares interpolation formula in the ver-
sion of Pohánka (2005) is the simplest possible function that
meets a set of basic requirements that an interpolation function
has to fulfil. The polynomial in local coordinates around the
point whose value it has to interpolate should be preferably of
the second order. The weights of the squared interpolation
errors at the measurement sites decrease with the distance to
the interpolation point. The behaviour of the function can be
controlled by two parameters: the smoothing distance d
0
and
the regularization distance d
1
.
From stress to strain colouring
The choice of the tensor field visualization method depends
on how much we want to show in order to be understood, as
well. Some methods require from the reader of the visualized
tensor field a sufficiently deep understanding of the tensor
notion and nature, or a superb 3D spatial imagination and
orientation, which can be a quite high expectation for some
users. A simple single-indexed colouring of stress states,
whose typical geophysical application is the WSM (Heidbach
et al. 2009), is on the other end of the stress state visualization
methods’ spectrum. Bada et al. (2007) goes one step further
and interpolates the discrete colours attributed to the stress
state points into a continuous stress state colour map. The less
colours the stress mapper has available, the more industrious
(s)he has to be in design of the glyphs for coding the stress
states. Still, a greyscale stress map remains a big challenge for
the reader (Jarosiński et al. 2006).
Our stress colouring method, described in Bednárik &
Kohút (2012) in detail, and now adapted to colouring strains
as well, aims to produce a reader-friendly, yet infor mation-
loaded colour map of a stress or strain field. The stress or
strain colour map however differs from the standard colour
map of a scalar quantity like height or temperature in one
important aspect: our colours depend on three scalar quantities
(the three principal stress or strain values) at once.
The stress tensor colouring method introduced in Bednárik
& Kohút (2012) needs an adjustment to the current applica-
tion. The main challenge is the unavailability of the relevant
rock characteristics (elastic and rheological parameters, mate-
rial strength, planes of weakness etc.) at the computation
points. In a complex geological environment to be expected in
the Tatra region, the stress field is influenced by many
variously oriented faults, some well exposed and documented,
some unknown. In this situation, the construction of a surface
stress field would be a too daring enterprise. Instead, we rather
opted for adaptation of the original method to the colouring of
the surface strain rate tensor field.
For the sake of derivation of the formulae, let us start with
the simplest isotropic elastic model. Here, the Hooke’s law in
Cartesian coordinates has the form
,
2
2
2
0
0
0
0
0
0
2
0
0
0
2
0
0
0
2
2
1
1
3
3
2
3
3
2
2
1
1
2
1
1
3
3
2
3
3
2
2
1
1
+
+
+
=
e
e
e
e
e
e
m
m
m
m
l
l
m
l
l
l
m
l
s
s
s
s
s
s
(5)
(
)
∑∑
=
=
=
M
m
N
n
n
m
n
m
I
y
x
c
y
x
v
0
0
,
513
SURFACE STRAIN RATE COLOUR MAP OF THE TATRA MOUNTAINS
GEOLOGICA CARPATHICA
, 2016, 67, 6, 509 – 524
where
l
and
m
are the Lamé elastic parameters at the point.
The symmetrical lower diagonal part of the stiffness matrix
has been omitted.
Our original colouring method was based on transforming
the principal stress values into so called colourables D, M, A or
their rescaled values d, m, a. We are interested to see whether
similar colouring can be done with the principal strains. The
first to be noted consequence of the isotropic Hooke’s law (5)
is that the principal stress directions and the principal strain
directions in isotropic elastic model are coaxial (we leave the
proof to the reader).
Let us have an isotropic elastic halfspace limited by a hori-
zontal free surface. Within the surface, let us introduce the
horizontal coordinate axes and denote them e
1
, e
2
. The remaining
coordinate axis e
3
will be vertical. If the horizontal surface is
traction-free, then
s
33
= 0,
s
23
= 0,
s
31
= 0, and the relationship of
the relevant stresses and strains degenerates to a 2D case
(called plane stress state), described by the matrix equation
(
) (
)
(
)
(
)
(
) (
)
+
+
+
+
+
+
=
2
1
2
2
1
1
2
1
2
2
1
1
2
0
0
0
2
4
2
2
0
2
2
2
4
e
e
e
m
m
l
m
l
m
m
l
m
l
m
l
m
l
m
l
m
l
m
s
s
s
,
(6)
or, after the substitution of
(
)(
)
ν
ν
ν
l
2
1
1
-
+
=
E
,
(
)
ν
m
+
=
1
2
E
,
where E is the Young’s modulus and
ν
is the Poisson’s
constant
(
)
-
-
=
2
1
2
2
1
1
2
2
1
2
2
1
1
2
2
1
0
0
0
1
0
1
1
e
e
e
ν
ν
ν
ν
s
s
s
E
.
(7)
It is obvious that two of the principal stress and strain direc-
tions are horizontal (generally not equal to e
1
, e
2
) and the third
is vertical i.e. equal to e
3
. Let us adopt the notation
s
1
,
s
2
(
e
1
,
e
2
) for horizontal principal stresses (strains) and
s
3
(
e
3
) for the
vertical principal stress (strain). Henceforth we completely
abandon the widely used convention of sorting and indexing
the principal stresses and strains according to their size. At the
horizontal traction-free surface, the principal stress value
s
3
corresponding to the vertical principal direction equals to zero.
This assumption is valid for any constitutive model.
In the isotropic elastic case, one of the principal directions
of strain is vertical, as well as with stress. Nevertheless, the
vertical principal strain value
e
3
is in general non-zero:
(
) (
)
(
) (
)
ν
e
e
ν
m
l
e
e
l
e
-
+
-
=
+
+
-
=
1
2
2
1
2
1
3
. (8)
After some algebra, we discover two useful relationships
between the principal stress and the principal strain values
valid in the case of plane stress:
,
(9)
(
)
(
)
(
) (
)
(
)
ν
e
e
e
e
m
s
s
t
+
-
=
-
=
-
=
1
2
2
2
1
2
1
2
1
E
m
. (10)
(
)
(
) (
) (
) (
) ( )
(
)
ν
e
e
m
l
m
l
m
e
e
s
s
s
-
+
=
+
+
+
=
+
=
1
2
2
2
3
2
2
1
2
1
2
1
E
m
The horizontal mean stress
s
m
and the horizontal maximum
shear stress
t
m
are the basic stress tensor shape descriptors in
the case of plane stress (or plane strain) or can be combined to
produce some more advanced shape descriptors — for instance
the descriptor (15).
The dependence of
e
3
on
e
1
,
e
2
(8) as well as the reducibility
of the plane stress state to a 2D horizontal problem (6), (7) are
very welcome facts to us. Having in mind the large error of the
vertical GNSS velocity component, confirmed also by Papčo
(2010) during his field experiments in the Tatra Mts., we can
now discard the vertical GNSS data and take into account only
the horizontal velocity components and compute the horizontal
2D strain or stress rate tensor. The eigenvalues
1
e
,
2
e
(
1
s
,
2
s
)
can then be completed by
3
e
(
3
s
) to the full 3D strain or stress
rate state.
Let us seek all colourable stress tensor descriptors as func-
tions
f
with the property
(
)
(
)
3
2
1
3
2
1
,
,
,
,
e
e
e
s
s
s
f
C
f
=
,
(11)
where C is a scaling factor. The scaling factor will then be
eliminated in the process of rescaling the output range of
f
to the interval 〈-1, 1〉. Thus, the procedures for compu ting such
tensor characteristics will stay the same and will yield exactly
the same colours, whether applied to principal stresses or
strains. Most importantly, the strain states can be coloured
without (almost) any knowledge of material parameters,
which is our goal.
We immediately realize that the descriptors
s
m
(9),
t
m
(10)
have exactly the desired form (11). The maximum shear stress
t
m
is the quantity upon which is based the Tresca criterion, not
the most appropriate to be used for the brittle surface of the
Earth’s crust. The problem with the more realistic Mohr-
Coulomb criterion
D
c
=
|
t
m
| - (
S
0
cos
φ
+
s
m
sin
φ
)
(12)
whether describing the failure on nascent or pre-existing
faults, is that it cannot be factorized like (11). Thus, the
Mohr-Coulomb (unlike the Tresca) criterion D
c
cannot be used
for colouring of strain states without the knowledge of the
relevant material parameters. We have to be satisfied with the
strain-compatible Tresca criterion as a rough indicator D of
the proneness to failure (a not very different alternative could
be the von Mises yield criterion).
In the plane stress case at the surface, the vertical principal
stress always equal to zero narrows the range of the values,
that the relevant mean stress
s
m
controlling the fault friction in
the Mohr-Coulomb criterion (12) can take, to the vicinity of
zero. With
s
m
not far from 0, the difference of the Tresca or
von Mises criterion and the Mohr-Coulomb criterion becomes
less dramatic and the Tresca (or von Mises) criterion becomes
an acceptable option.
The last step in computing the first colourable d is the nor-
malization of D∈〈D
min
, D
max
〉
to d∈〈-1, 1〉:
514
BEDNÁRIK, PAPČO, POHÁNKA, BEZÁK, KOHÚT and BRIMICH
GEOLOGICA CARPATHICA
, 2016, 67, 6, 509 – 524
≠
-
-
-
=
=
min
max
min
max
min
max
min
max
for
2
for
0
D
D
D
D
D
D
D
D
D
d
,
, (13)
where
(
)
3
2
3
1
2
1
,
,
max
e
e
e
e
e
e
m
-
-
-
=
D
(Tresca failure
criterion used), and D
max
, D
min
are the maximum and minimum
of D to be shown.
The best candidate for the colourable M, indicating whether
the stress/strain state is compressive or tensile, is the mean of
the horizontal principal stresses/strains
s
m
(9). The normaliza-
tion of M∈〈M
min
, M
max
〉
to m∈〈-1, 1〉 is analogical to (13):
≠
-
-
-
=
±
=
min
max
min
max
min
max
min
max
for
2
for
0
M
M
M
M
M
M
M
M
M
m
,
(14)
where M =
s
m
, and M
max
, M
min
are the maximum and mini-
mum of M to be shown. Like in (13), M
max
, M
min
do not need
to be the actual extrema of M within the investigated region. In
that case, however, colours generated for M > M
max
or M < M
min
are the same as the colours for M = M
max
or M = M
min
, respec-
tively — the palette used in the map becomes saturated.
Unlike in (13), it is convenient to map M = 0 to m = 0.
To assure this, M
min
has to be set to M
min
= – M
max
(or M
max
to
M
max
= – M
min
). With the ± switch we added the possibility to
invert the scale.
The last colourable A (rescaled to a) shall be a predictor of
the fault orientation. In its original form (Bednárik & Kohút
2012), a has not made a clear colour distinction between the
regions with different orientations of the theoretical fault
plane. Here we propose the solutions. Let us show the fault
planes, whose normal vector has the vertical component equal
to zero by cold colours and all other fault planes by warm
colours (or vice versa). If
s
1
<
s
3
<
s
2
or
s
2
<
s
3
<
s
1
, the stress
state would favour the horizontal strike-slips on fault planes,
whose normal vectors are horizontal. Here once again, we take
advantage of
s
3
= 0. Then,
|
t
m
|
> |
s
m
|
means that the horizontal
principal stresses
s
1
,
s
2
have different signs and
s
3
= 0 falls
between them. Thus,
(
)
(
)
(
)
(
)
ν
e
e
ν
e
e
m
l
e
e
m
l
e
e
s
t
+
+
-
-
=
+
+
+
-
=
=
1
1
2
3
2
2
1
2
1
2
1
2
1
m
m
A
(15)
is a suitable indicator of the fault orientation, because it
also meets the condition (11). One complication arises from
the need to map the output range of A to 〈-1, 1〉, the other from
the fact that a fault orientation indicator has to be scaled
locally, depending only on the strain eigenvalues at the point.
The scaling factor (1– v) /(1+ v) cannot be eliminated by the
norma lization over the whole region, unlike in the case of D
or M. The only way to eliminate the scaling factor is to esti-
mate it based on the estimate of the Poisson’s constant. Fortu-
nately, its variations in the relevant crustal rocks are not large
in relation to the goals of this study. The value of the scaling
factor can be guessed and eliminated with sufficient
accuracy.
An elegant way to map A∈〈0, ∞) in (15) to A
log
∈(-∞
, ∞) and
A = 1 to A
log
= 0
is to take the logarithm (of any suitable base b)
(
)
(
)
+
-
+
+
-
±
=
±
=
ν
ν
e
e
e
e
s
t
1
1
log
log
log
2
1
2
1
log
b
b
m
m
b
A
(16)
and normalize subsequently to a∈〈-1, 1〉 (similarly to (14)),
with the option to cut off the extreme values. In real data, we
almost never hit the singularities of A or A
log
, nevertheless, we
have to prevent the overflows. With the
±
switch we added the
possibility to invert the scale, just as in the case of m.
We can offer a different and more general solution to the
fault plane descriptor of the class (11). The term
(
)
2
3
3
1
2
1
s
s
s
s
s
s
-
+
-
-
(17)
is singular for the hydrostatic stress state
s
1
=
s
2
=
s
3
. In all
other cases it evaluates to 1, if the vertical principal stress
value
s
3
falls between the horizontal eigenvalues
s
1
and
s
2
(or
is equal to either of the two), otherwise yields a value from the
range 〈0, 1). Since we want to have a colour gradation in the
first case, as well, and want to cover the whole output range
a∈〈-1, 1〉, we define a as a forked function
≠
-
+
-
-
-
-
+
-
-
=
-
+
-
-
-
-
±
=
1
for
1
1
for
2
3
3
1
2
1
2
3
3
1
2
1
2
3
3
1
2
1
2
1
2
3
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
s
a
.
.
(18)
The term (
s
3
–
s
2
) / (
s
1
–
s
2
) in the upper branch is formally
identical with the stress shape ratio
F
used in structural geo-
logy. The descriptor (18) can be used not only in the case of
plane stress on the free surface, as (15), (16), but also in a more
general case of
s
3
≠
0, but still considerable as vertical prin-
cipal stress. In the most general application, from the three
principal stress directions, one could be selected as “the most
vertical” and the respective principal stress values can be
denoted accordingly to comply with the convention behind
(18). In favour of the convertibility to the strain formulation, it
is also advantageous not to substitute
s
3
= 0 even in the case of
plane stress. The relationship (
s
i
–
s
j
) / 2=
m
(
e
i
–
e
j
) is valid for
any mutually corresponding couples of stress and strain prin-
cipal values
s
i
,
e
i
and
s
j
,
e
j
(and for i = j, as well). Then in
terms of strains, (19) is as simple as the original formula in
stresses (18):
≠
-
+
-
-
-
-
+
-
-
=
-
+
-
-
-
-
±
=
1
for
1
1
for
2
3
3
1
2
1
2
3
3
1
2
1
2
3
3
1
2
1
2
1
2
3
e
e
e
e
e
e
e
e
e
e
e
e
e
e
e
e
e
e
e
e
e
e
a
.
(19)
515
SURFACE STRAIN RATE COLOUR MAP OF THE TATRA MOUNTAINS
GEOLOGICA CARPATHICA
, 2016, 67, 6, 509 – 524
The only problem is that in plane stress,
e
3
≠
0. Thus, the
e
3
given by (8) and involving the Poisson’s constant has to be
substituted (we leave the task to the reader). Even in this case,
we ended up with the necessity to enter a single material
parameter, which we wanted to avoid.
At this point, we have two remarks. In A=
|
t
m
|
/ |
s
m
|
, we
have an obvious example that out of the three tensor descrip-
tors D, M, A of the plane stress state, only the first two are
independent. This corresponds to the degenerate 2D character
of the plane stress state, where only the horizontal principal
strain values are independent and the vertical one is their
linear combination. In spite of this, the colourable a plays an
important role in the colour generation. A necessary conse-
quence of the dependence of A on D and M is a somewhat
reduced colour space compared to the full colour cube
(Fig. 3a, b), which is only a minor flaw. In a according to (19),
that interdependence of colourables is hidden, but still present.
The other remark applies to the now unveiled fact that
expressing the stress tensor descriptors in terms of principal
strains is a quite unnatural task. The more the descriptor
pertains to the balance of forces — like in the case of Mohr-
Coulomb failure criterion or the fault plane preference indicator
— the bigger the need for the input of material parameters to
convert the strains back to stresses. We tried to keep that need
as small as possible, but in the case of A, we had to make an
exception and admit the Poisson’s constant to enter our
calculations.
The herewith introduced tensor shape characteristics D, M,
A and their rescaled values d, m, a can be regarded in two
different ways: we either take them for stress tensor shape
descriptors expressed in terms of principal strains, under the
assumption of isotropic elasticity (or some other simple con-
stitutive model); or we completely ignore their stress rootage
and any underlying assumptions and take them for indepen-
dent strain state characteristics. They can then be used as strain
rate tensor descriptors, as well.
To generate the colours corresponding to the colourables d,
m, a, the user can choose from several RGB colour function
formulae we constructed:
(
)
(
)
(
)
T
d
a
d
a
d
m
a
m
d
d
a
m
d
a
d
m
a
m
d
d
a
m
d
a
a
m
d
c
-
-
-
+
-
-
+
-
-
+
+
-
+
+
+
-
+
=
3
1
3
1
3
4
1
,
,
(20)
generates stress colours on white background, where white
(1, 1, 1) is associated with minimum Tresca stress (d = –1).
(
)
(
)
(
)
T
Blk
d
a
d
a
d
m
a
m
d
d
a
m
d
a
d
m
a
m
d
d
a
m
d
a
a
m
d
c
-
+
-
+
-
-
+
-
+
+
+
-
+
+
+
+
+
=
1
1
1
1
1
4
1
,
,
(21)
generates stress colours on black background, where black
(0, 0, 0) is associated with minimum Tresca stress (d = –1).
The simpler, however less vivid colours generating trilinear
stress colour functions are also available for both
backgrounds:
(
)
T
lin
d
a
d
a
m
d
m
d
m
d
m
d
a
m
d
c
-
-
-
-
-
-
+
+
-
=
3
3
3
4
1
,
,
,
(22)
(
)
T
Blk
lin
d
a
d
a
m
d
m
d
m
d
m
d
a
m
d
c
-
+
-
-
-
+
+
+
+
=
1
1
1
4
1
,
,
.
(23)
With these settings, we produce a three-dimensional colour
palette changing from green to red as the strain/stress state
changes from extensional to compressional, from dark (palette
(21)) or pale (palette (20)) to bright colours as the strain/stress
state gets closer to failure, from warm to cold hues as the
theoretical fault planes become more vertical and more strike-
slip. In our opinion, this colour convention is quite logical
and intuitive.
With a change of sign in (14), (16), (19) we could have pro-
duced a palette more similar to the convention of the World
Stress Map (where normal fault is red, strike-slip green and
thrust fault blue), but we deliberately chose not to do so. One
reason is that we wanted to underline the novelty and inde-
pendence of our colouring method. The other reason is that
due to the logic of the hue – saturation – value (HSV) colour
system, we are not able to adhere to WSM convention totally,
anyway. To avoid the risk of misleading the reader used to
the WSM, we rather opted for a completely different colouring
scheme.
Colour maps of strain rates
The colour maps presented here need some translation into
human speech in order to be fully appreciated. We are able to
distinguish by colours the basic tectonic regime categories:
1. warm red — shortening — thrusting;
2. cold red — transpression;
3. blue — vertical strike-slip faulting;
4. cold green — transtension;
5. warm green — extension — normal faulting or vertical
tearing.
In the colour maps, these categories should be understood as
likely tendencies, not as sharp statements.
Another help for the reader’s imagination can be found in
the three-dimensional colour palettes corresponding to (20)
and (21) (Fig. 3a, b). They will be used for colouring all the
maps in this chapter. In the captions to the figures with maps,
we will only specify the respective ranges of D, M, A that
were normalized to d∈〈-1, 1〉, m∈〈-1, 1〉, a∈〈-1, 1〉.
Once we have grasped the basic logic of the stress colou-
ring, we can start to enjoy the reading of combined map of
stress, faults and elevations. The possibility of a relatively
easy multi-contextual map reading is one of the main advan-
tages of our stress colouring scheme — in the World Stress
Map, already the stress layer alone is hardly readable.
516
BEDNÁRIK, PAPČO, POHÁNKA, BEZÁK, KOHÚT and BRIMICH
GEOLOGICA CARPATHICA
, 2016, 67, 6, 509 – 524
Let us start with a map for which we do not need any colours
at all (Fig. 4). The measured velocities are fitted by a purely
translational field, that generates zero strain rates. In the case
of the whole-region least squares fit, the maximum approxi-
mation error was Dv
y
= – 0.712 mm/year in LOMS (greater
than the respective standard deviation). In the case of the
whole-region adaptive weights fit, the maximum approxi-
mation error was Dv
y
= 0.602 mm/year in ZISE (greater than
the respective standard deviation). In Fig. 4, we show only the
more successful adaptive weights fit.
a
b
Fig. 3. a — 3D colour palette on white background (20) used for colouring the stress/strain (rate) states in front view (left) and rear view (right).
b — 3D colour palette on black background (21) used for colouring the stress/strain (rate) states in front view (left) and rear view (right).
517
SURFACE STRAIN RATE COLOUR MAP OF THE TATRA MOUNTAINS
GEOLOGICA CARPATHICA
, 2016, 67, 6, 509 – 524
Let us continue with the simplest possible model that gene-
rates nonzero strain rates and thus, colours different from the
background. It is the bilinear velocity field ((2) for M = 1,
N = 1). In the case of the whole-region least squares fit, the
maximum approximation error was Dv
y
= – 0.560 mm/year in
LOMS (greater than the standard deviation). In the case of the
whole-region adaptive weights fit, the maximum approxima-
tion error was Dv
y
= 0.487 mm/year in TAKO (slightly smaller
than the standard deviation).
We see that the whole-region bilinear velocity field can
model considerable variations of the velocity directions and in
spite of that, generate a strain field, whose tensor descriptors
vary only moderately. Moreover, the velocity fields that both
are acceptable as the solutions of the interpolation problem
generate contradictory colour maps (Fig. 5 a, b).
Our preliminary conclusion from this contradiction is that
on a regional scale, no relevant strains occur within the speci-
fied region of the Tatra Mts. The search for relevant strains
should be focused on more detailed scales and local structures.
To sharpen the lens, we have to increase the polynomial order
in (2). The problem with this approach is that discarding
GANO and LIES because of their unsatisfactory geological
setting, we are left with 9 data points. The number of unknown
coefficients of the biquadratic fit is exactly 9. The system is
just-determined and has a unique solution, that is however no
reason for jubilation. In this case, we no longer have a minimi-
zation problem, but in its stead an exact interpolation, that
strictly adheres to the possibly unreliable values at the data
points. The interpolated velocity field has a huge ratio of the
largest to the smallest velocity and the velocity directions
behave rather unbelievably
between the data points that are
close neighbours (Fig. 6a).
One — mathematically not very
satisfactory — remedy would be
an incomplete polynomial, that
grants the fit some freedom. If we
discard the highest powers term
c
22
x
2
y
2
in the biquadratic polyno-
mial, we get a smoother velocity
and strain field (Fig. 6b). In both
cases, the interpolation errors at
the measurement points are so
small that it makes no sense to
start the iterations of the adaptive
weights method.
Increasing the polynomial order
in the whole-region approxi-
mation allows it, allegorically
speaking, to develop “weird fanta-
sies” on why the two data points
next to each other have different
velocities, because it is misad-
vised by very distant data points
that have a very distant idea about
the true relationship of the former
two data points.
Therefore, the distant data points should not be taken as
seriously as the data points next to the interpolation point. This
is the main idea behind the moving least squares method.
Because of the locally valid fit of the data, we label it as a local
interpolation, in contrast to the whole-region interpolation. As
a matter of fact, the moving least squares method in the imple-
mentation of Pohánka (2005) takes into account the complete
set of data points just as the whole-region methods, but for
each interpolation point, it generates its own interpolation for-
mula, valid only at that single point.
The results of the moving least squares (Fig. 7b) are in good
accordance with the biquadratic whole region fits (Fig. 6 a, b).
Obviously, only after going into the second polynomial order
and into the local scales, the interpolation methods are able to
focus on the Sub-Tatra fault system and parallel structures,
whose motions were averaged out in the bilinear whole
region fit.
For the reasons of the dependence of the colourable a on
d and m and also because of the palette saturated in those
maps, where we had to ignore the extreme strain artefacts
at the margins of the region, the colours generated by our
method resemble pure spectral colours. This is an advantage
that can be used in the diagnostics of differences between
two (or more) maps: where very different colours are super-
imposed onto each other, the resulting blend is a dirty,
khaki colour that can be easily recognized in the com pound
map. To assure the fairness of the comparison, all the
palette range settings in all thus compared maps have to be
the same.
Fig. 4. Measured (black arrows) and interpolated (grey arrows) surface velocities shown together
with height contours (interval 200 m) and fault scheme according to Bezák et al. (2008). The two
velocity fields have different scales. Adaptive weights least squares interpolation of the 0-th order
was used. The stations included in the interpolation are represented by white dots, the excluded
stations by black dots.
518
BEDNÁRIK, PAPČO, POHÁNKA, BEZÁK, KOHÚT and BRIMICH
GEOLOGICA CARPATHICA
, 2016, 67, 6, 509 – 524
b
a
Fig. 5. Strain rate colour map, measured (yellow arrows) and interpolated (grey arrows) surface velocities shown together with height contours
(in orange, interval 200 m) and fault scheme according to Bezák et al. (2008). The two velocity fields have different scales.
a — Whole- region least-squares bilinear interpolation was used. The vertices of the black background palette ((21), Fig. 3b) correspond to
D
min
= 15 nstrain/year, D
max
= 20 nstrain/year, M
min
= – 5 nstrain/year, M
max
= 5 nstrain/year,
log
min
A
= – 1.74,
log
max
A
= 1.74 (logarithmic scale (16) with
b = 2 used). A part of the D, M, A triplets comes out of the specified range and is projected onto the surface of the palette cube (palette
saturation).
b — Whole- region adaptive weights least-squares bilinear interpolation was used. The vertices of the black background palette ((21), Fig. 3b)
correspond to D
min
= 0 nstrain/year, D
max
= 20 nstrain/year, M
min
= – 20 nstrain/year, M
max
= 20 nstrain/year,
log
min
A
= – 1.27,
log
max
A
= 1.27 (logarithmic
scale (16) with b = 2 used). A part of the D, M, A triplets comes out of the specified range and is projected onto the surface of the palette cube
(palette saturation).
519
SURFACE STRAIN RATE COLOUR MAP OF THE TATRA MOUNTAINS
GEOLOGICA CARPATHICA
, 2016, 67, 6, 509 – 524
b
a
Fig. 6. Strain rate colour map, measured surface velocities (yellow arrows) and interpolated surface velocity directions (grey arrows) shown
together with height contours (in orange, interval 200 m) and fault scheme according to Bezák et al. (2008).
a — Whole-region least-squares fully biquadratic interpolation was used. The vertices of the black background palette ((21), Fig. 3b) corre-
spond to D
min
= 0 nstrain/year, D
max
= 100 nstrain/year, M
min
= – 100 nstrain/year, M
max
= 100 nstrain/year,
log
min
A
= – 1.22,
log
max
A
= 1.22 (logarithmic
scale (16) with b = 2 used). A considerable part of the D, M, A triplets (mainly at the region’s margins) comes out of the specified range and is
projected onto the surface of the palette cube (palette saturation).
b — Whole-region least-squares incomplete biquadratic interpolation was used. The vertices of the black background palette ((21), Fig. 3b)
correspond to D
min
= 0 nstrain/year, D
max
= 100 nstrain/year, M
min
= – 100 nstrain/year, M
max
= 100 nstrain/year,
log
min
A
= – 0.422,
log
max
A
= 0.422 (loga-
rithmic scale (16) with b = 2 used). A considerable part of the D, M, A triplets (mainly at the region’s margins) comes out of the specified range
and is projected onto the surface of the palette cube (palette saturation).
520
BEDNÁRIK, PAPČO, POHÁNKA, BEZÁK, KOHÚT and BRIMICH
GEOLOGICA CARPATHICA
, 2016, 67, 6, 509 – 524
b
a
Fig. 7. Strain rate colour map, measured surface velocities (yellow arrows) and interpolated surface velocity directions (grey arrows) shown
together with height contours (in orange, interval 200 m) and fault scheme according to Bezák et al. (2008).
a — Moving least-squares interpolation (smoothing distance d
0
= 10 km, regularization distance d
1
= 320 km) was used with all 11 stations
involved. The vertices of the white background palette ((20), Fig. 3a) correspond to D
min
= 0 nstrain/year, D
max
= 100 nstrain/year,
M
min
= – 100 nstrain/year, M
max
= 100 nstrain/year,
log
min
A
= – 1.0,
log
max
A
= 1.0 (logarithmic scale (16) with b = 2 used). A considerable part of the D,
M, A triplets (mainly at the region’s margins) comes out of the specified range and is projected onto the surface of the palette cube (palette
saturation).
b — Moving least-squares interpolation (smoothing distance d
0
= 10 km, regulari zation distance d
1
= 320 km) was used with 9 stations involved.
The vertices of the white background palette ((20), Fig. 3a) correspond to D
min
= 0 nstrain/year, D
max
= 100 nstrain/year, M
min
= – 100 nstrain/year,
M
max
= 100 nstrain/year,
log
min
A
= – 1.0,
log
max
A
= 1.0 (logarithmic scale (16) with b = 2 used). A considerable part of the D, M, A triplets (mainly at the
region’s margins) comes out of the specified range and is projected onto the surface of the palette cube (palette saturation).
521
SURFACE STRAIN RATE COLOUR MAP OF THE TATRA MOUNTAINS
GEOLOGICA CARPATHICA
, 2016, 67, 6, 509 – 524
We demonstrate it on the example of two maps generated
by the moving least squares interpolation. In the first one, the
interpolation is based on all 11 original data points, including
GANO and LIES. In the second one, the interpolation is
based on 9 data points without GANO and LIES. We clearly
see that the khaki spots cover not only the areas next to
GANO and LIES, but all the areas insufficiently covered by
data points.
Discussion
The detection of the Sub-Tatra fault system in the second
order interpolation by such a low number of GNSS stations
(even not deployed across the fault) is pleasant and suspicious
at the same time. It may be argued that the line separating the
compression and extension zone is an artefact determined by
the geometry of the GNSS network rather than by the true geo-
dynamics. Another problem is the reliability of the biquadratic
approximation that does not (because it cannot) take into
account the high data uncertainties. The interpolated velocity
field in the band between the two groups of stations does not
behave according to the natural expectations. The decision not
to pay attention to the changing principal directions of strain
rates and display only the tensor shape descriptors was correct,
given the data uncertainties and interpolation errors visible in
the interpolated velocity field.
On the other hand, we have independent evidence consistent
with the extension of the southern part and compression of the
northern part of the Tatra Mts. that we have found in this study.
It is a known fact that the neotectonic uplift of the Tatra horst
was asymmetrical, with faster uplift rates in the south and tilting
down towards the north. This is the reason why the highest
Tatra peaks are found southwards from the main ridge.
Recently, Králiková et al. (2014) studied the tectonic evolu-
tion of the Tatra Mts. based on combined structural, sedimen-
tary, geomorphological, and fission track data evidence.
The study by Anczkiewicz et al. (2015) confirmed the asym-
metry of the Tatra Mts. uplift by low-temperature thermochro-
nology. The study by Szalaiová et al. (2008) provides
a gravimetrical indication to the southward dipping of the
Sub-Tatra fault.
Whereas we take for trustworthy the compression in the
north and extension in the south of the Tatra Mts. computed
from the GNSS data, the line between the two zones deserves
some further discussion. It is parallel to the Sub-Tatra fault
system. Nevertheless, if the extension and compression should
be linked to the tilting of the horst, it should be rather identi-
fied with the axis of the tilt, and the Sub-Tatra fault system
rather with the warm green maxima. The exact course of the
two lines cannot be determined from that sparse and uncertain
GNSS data.
The Sub-Tatra fault could have eased the post-glacial
rebound of Tatra Mts. In this case, the strain rate pattern along
the normal fault should be extensional, as well. If the faults on
the northern foot of the Tatra Mts. are harder to mobilize than
the Sub-Tatra fault, the uplift due to the rebound could be
tilted northwards, just as the neotectonic uplift.
The orders 0, 1 and 2 of the interpolating polynomials can
be interpreted as levels of trust in the surface velocity GNSS
Fig. 7. c — Strain rate colour blend map and measured surface velocities (yellow arrows) shown together with height contours (in orange,
interval 200 m) and fault scheme according to Bezák et al. (2008). The colour map for 11 stations (Fig. 7a) and 9 stations (Fig. 7b) were blended
together. Where the maps are not in mutual agreement, colours different from the colours in (Fig. 7a, b) appear. The bigger the disagreement,
the dirtier the resulting blend.
c
522
BEDNÁRIK, PAPČO, POHÁNKA, BEZÁK, KOHÚT and BRIMICH
GEOLOGICA CARPATHICA
, 2016, 67, 6, 509 – 524
data and, at the same time, as degrees of cartographic gene-
ralization coupled to the scale of observation.
The order 0 corresponds to the deep mistrust that us allows
to interpret the GNSS data only as pure translation (although
this mistrust is slightly deeper than the standard deviations of
the data permit). In the resulting map, no strain/stress rates
appear. In the scales of whole-world stress/strain mapping, the
stress/strain rates in the region are really negligible.
The order 1 corresponds to a cautious and balanced approach
to the velocities and their deviations. This order allows for
considerable variations of the data within their uncertainties.
The two colour maps generated from the interpolated data
contradict each other. The order 1 also corresponds to the
regional scale cartographic generalization, in which the local
behaviour is averaged out in the regional trends. These proved
insignificant, as well.
The order 2 combined with the low number of stations
mani fests an undue trust in the velocity means and undue
ignorance of their uncertainties. On the other hand, only this
and higher orders allow to focus on the local scale behaviour
that got oversmoothened in the order 1 interpolation. In our
application, however, the order 2 proved unable to preserve
realistic regional trends alongside the local behaviour.
One viable solution we found for going into local scales was
a local interpolation method – the moving least squares
(Pohánka, 2005). With only 9 stations that qualified into the
selection, the performance of this method was far from the
optimum that would be achievable in an application with more
sites more evenly distributed.
Thus, all the solutions provided by the interpolations of
order 0, 1, 2 (or rather the moving least squares in the latter
case) are possible answers to the question of what the GNSS
data tell us about the surface deformations in the Tatra Mts.
region. The reader’s choice of the right answer depends on his
trust in the data and on the spatial scale he is interested in.
The colour rendering of the strain rates allows an easy visual
comparison of the strain rate maps generated for different
parameter settings or different input data. A simple average of
the RGB values in the maps efficiently displays their diffe-
rences as dirty spots. The compound map would thus be
automatically partitioned into subregions with different
degrees of reliability. In Mentlík & Novotná (2010), the scien-
tific reliability rating of the map content is based on an expert’s
assessment. Our two maps prepared by moving least squares
interpolations for 11 and 9 stations are in mutual accordance
mainly in the area with higher density of the stations. The rest
of the region shows an instability of the interpolation solution
(although exaggerated by removal of the two most proble-
matic stations GANO and LIES). Ideally, the stations should
be sufficiently numerous to ensure that removing a few of
them does not change the colour map very substantially.
The colour maps averaging can be used in an effective
visuali zation of the data uncertainties, as well. For instance,
we could repeatedly generate random site velocities with the
same mean and same standard deviation as in the measured
dataset. For each randomly generated dataset, we could
compute the colour map and finally collect the maps, average
them and partition them into areas of different reliability.
Alternatively, the time span of the original measurements
could be split into subintervals and for each of them, a dataset
of stations’ short term average velocities could be computed.
The evolution of the strain field and of its uncertainties in time
could be animated by a sequence of individual colour maps or
by adding map by map to the average compound colour map.
Such an animation is, however, difficult to show in the form of
a classical journal paper.
Conclusion
The horizontal surface velocities and their standard devia-
tions at 9 or 11 stations within the Tatra Mts. region allow us
to construct a broad range of interpolated velocity field solu-
tions. The colour maps of the whole-region least squares
approximations up to the polynomial order 1 show no remar-
kable regional trends of surface strain rates. The whole-region
least squares fits by higher order polynomials are problematic
not only because of the low number of data points. Thus, to
resolve more local variations of the strain rates, local interpo-
lation methods should be used. Out of them, we tested the
moving least squares method and it yielded promising results.
For the full deployment of the local interpolation methods’
potential, more (and at the same time, more evenly distri buted)
local data is needed: the most welcome would be the velocities
of site couples across the known faults, which are not abun-
dant in the current dataset.
The tensor colouring method now adapted to the strain rate
field at the Earth’s surface proved to be a useful tool in showing
the spatial variations of the relevant strain rate tensor shape
descriptors in the form of colour maps. The application of the
method even to a sparse dataset has clearly demonstrated its
advantages: modest computational requirements, high reso-
lution, easy intuitive readability of resulting maps, natural,
easy to implement and reader-friendly colour imaging of data
uncertainties, integrability into multi-contextual graphic
rendering of geophysical data by standard geographical infor-
mation systems.
Acknowledgements: The authors are grateful to the Slovak
grant agency VEGA (grants No. 2/0002/12 , 1/0642/13,
1/0141/15, 1/0714/15, 2/0042/15, 2/0091/15) and to the
Slovak Research and Development Agency (grants
APVV-0724-11 and APVV-0212-12) for the partial support
of this work.
References
Altamimi Z., Collilieux X., Legrand J., Garayt B. & Boucher C. 2007:
ITRF2005: New release of the International Terrestrial Refe-
rence Frame based on time series of station positions and Earth
Orientation Parameters. J. Geophys. Res. 112, B09401.
Anczkiewicz A.A., Danišík M. & Środoń J. 2015: Multiple
523
SURFACE STRAIN RATE COLOUR MAP OF THE TATRA MOUNTAINS
GEOLOGICA CARPATHICA
, 2016, 67, 6, 509 – 524
low-temperature thermochronology constraints on exhumation
of the Tatra Mountains: New implication for the complex evolu-
tion of the Western Carpathians in the Cenozoic. Tectonics 34,
11, 2296–2317.
Bada G., Horváth F., Dövényi P., Szafián P., Windhoffer G. &
Cloetingh S. 2007: Present-day stress field and tectonic inver-
sion in the Pannonian basin. Global Planet. Change 58, 1–4:
165–180.
Bednárik M. & Kohút I. 2012: Three-dimensional colour functions
for stress state visualization. Computers & Geosciences 48,
117–125.
Bezák V., Broska I., Ivanička J., Ivanička J., Polák M., Potfaj M.,
Buček J., Janočko J., Kaličiak M., Konečný V., Šimon L., Elečko
M., Fordinál K., Nagy A., Maglay J. & Pristaš J. 2008: New
edition of general geological map of the Slovak Republic in the
scale 1:200,000 [Prehľadná geologická mapa Slovenskej repub-
liky 1:200 000]. Ministry of environment of the Slovak Republic,
State Geological Institute of Dionýz Štúr, Bratislava.
Bezák V., Broska I., Ivanička J., Reichwalder P., Vozár J., Polák M.,
Havrila M., Mello J., Biely A., Plašienka D., Kaličiak M., Žec
B., Elečko M., Janočko J., Preszlényi M., Marko M., Maglay J.
& Pristaš J. 2004: Tectonic map of the Slovak Republic at
1:500,000 scale [Tektonická mapa Slovenskej republiky
1:500 000]. Ministry of environment of the Slovak Republic,
Bratislava. ISBN 80-88974-62-3 (in Slovak).
Briestenský M. & Stemberk J. 2008: Monitoring mikropohybov v
jaskyniach západného Slovenska. Acta Carsologica Slovaca, 46,
2, 77–83.
Briestenský M. 2008: Indicators of fault activity in the Brezovské part
of Malé Karpaty Mts. [Indikátory zlomovej aktivity brezovskej
časti Malých Karpát]. Geomorphologia Slovaca et Bohemica,
8, 1, 16–25 (in Slovak).
Burchart J. 1972: Fission-track age determinations of accessory apa-
tite from the Tatra Mountains, Poland. Earth Planet. Sci. Lett.
15, 4.
Czarnecki K. & Mojzeš M. 2000: Geodynamics of Tatra Mountains
— present status and development of the project. Report on
geodesy 7, Warszawa,123–124.
Dach R., Hugentobler U., Fridez P. & Meindl M. 2007: Bernese GPS
Software version 5.0, Astronomical Institute, University of Bern.
Drewes H. 2009: The actual plate kinematic and crustal deformation
model APKIM2005 as basis for a non-rotating ITRF.
In: Geodetic reference frames. Springer, Berlin, 95–99.
EPN 2016: EUREF permanent network. Available online at:
http://www.epncb.oma.be.
Fejes I. 2002: Consortium for Central European GPS Geodynamic
Reference Network (CEGRN): Concept, Objectives and Organi-
zation. Reports on Geodesy 61, 1, 15–21.
Fojtíková L., Vavryčuk V., Cipciar A. & Madarás J. 2010: Focal
mechanisms of micro-earthquakes in the Dobrá Voda seismo-
active area in the Malé Karpaty Mts. (Little Carpathians),
Slovakia. T
ectonophysics 492, 1–4, 213–229.
Hefty J. & Gerhátová Ľ. 2006: Site velocities from long-term epoch
GPS observations — Case study: Central Europe regional geo-
dynamics project 1994–2005. Acta Geodyn. Geomater. 3, 3
(143), 7–17.
Hefty J. 2004: Global positioning system in 4D geodesy [Globálny
polohový systém v štvorrozmernej geodézií]. STU, Bratislava.
ISBN 80-227-2027-5, 1–112 (in Slovak).
Heidbach O., Tingay M., Barth A., Reinecker J., Kurfeß D. & Müller
B. 2009: The World Stress Map based on the database release
2008, equatorial scale 1:46,000,000. Commission for the Geo-
logical Map of the World, Paris, doi: 10.1594/GFZ.WSM.
Map2009.
Hradilek L. 1984: Alpine surveying, trigonometric leveling and
3-dimensional terrestrial triangulation [Vysokohorská geodézie
— trigonometrická nivelace a trojrozměrná terestrická trian-
gulace]. Academia, Praha, 1–230 (in Czech).
IGS 2016: International GNSS Service. Available online at:
http://igscb.jpl.nasa.gov/.
Jarosiński M., Beekman F., Bada G. & Cloetingh S. 2006: Redistribu-
tion of recent collision push and ridge push in Central Europe:
insights from FEM modeling. Geophys. J. Int. 167, 2, 860–880.
Košťák B. 1969: A new device for in-situ movement detection and
measurement. Experimental Mechanics 9, 8, 374–379.
Králiková S., Vojtko R., Sliva Ľ., Minár J., Fügenschuh B., Kováč M.
& Hók J. 2014. Cretaceous–Quaternary tectonic evolution of the
Tatra Mts. (Western Carpathians): constraints from structural,
sedimentary, geomorphic, and fission track data. Geol. Carpath.
65, 4, 307–326.
Littva J., Hók J. & Bella P. 2015: Cavitonics: Using caves in active
tectonic studies (Western Carpathians, case study). J. Struct.
Geol. 80, 47–56.
Mentes Gy. 2008: Observation of recent tectonic movements by
extensometers in the Pannonian Basin. J. Geodynamics 45, 4–5,
169–177.
Mentlík P. & Novotná M. 2010: Elementary forms and “scientific
reliability” as an innovative approach to geomorphological map-
ping. Journal of Maps 6, 1, 564–583.
Nemčok J., Bezák V., Biely A., Gorek A., Gross P., Halouzka R.,
Janák M., Kahan Š., Kotański Z., Lefeld J., Mello J., Reichwalder
P., Rackowski W., Roniewicz P., Ryka W., Wieczorek J. &
Zelman J. 1994: Geological map of the Tatra Mts. at 1:50,000
scale [Geologická mapa Tatier 1:50 000]. Ministry of environ-
ment of the Slovak Republic , State Geological Institute of
Dionýz Štúr, MOSZNL, PIG, Bratislava (in Slovak).
Nemčok J., Bezák V., Janák M., Kahan Š., Ryka W., Kohút M.,
Lehotský I., Wieczorek J., Zelman J., Mello J., Halouzka R.,
Rackowski W. & Reichwalder P. 1993: Explanatory notes to the
geological map of the Tatra Mts. at 1:50,000 scale [Vysvetlivky
ku geologickej mape Tatier 1: 50 000]. State Geological Institute
of Dionýz Štúr, Bratislava, 1–135 (in Slovak).
Papčo J. 2010: Monitoring of the Earth’s Crust Deformations in
Alpine Area [Monitorovanie deformácií zemskej kôry vo vyso-
kohorskom prostredí]. Thesis, Slovak University of Technology
in Bratislava, Faculty of Civil Engineering, 1–127 (in Slovak,
with English summary).
Pohánka V. 2005: Universal interpolation method for many-dimen-
sional spaces. Available online at: http://gpi.savba.sk/GPIweb/
ogg/pohanka/int.pdf.
Sperner B. & Ratschbacher L. 1995: Rise and fall of the High Tatra
Mts. In: Proceedings of Europrobe Workshop Pancardi. SAV,
Bratislava, 67–69.
Stangl G. 2007: Guidelines for CEGRN Reprocessing version 1.3.
Available online at: http://cergops2.iwf.oeaw.ac.at.
Szalaiová E., Bielik M., Makarenko I., Legostaeva O., Hók J.,
Starostenko V., Šujan M. & Šefara J. 2008: Calculation of
a stripped gravity map with a high degree ofaccuracy: a case
study of Liptovská Kotlina Basin (Northern Slovakia). Geol.
Quarterly 52, 2, 103–114.
524
BEDNÁRIK, PAPČO, POHÁNKA, BEZÁK, KOHÚT and BRIMICH
GEOLOGICA CARPATHICA
, 2016, 67, 6, 509 – 524
Appendix
Surface deformation velocity in Tatra Mountains region — source data by Papčo (2010)
# code
X_east_UTM Y_north_UTM
Z DY_north
SDY
DX_east
SDX DZ_vertical SDZ_vertical
[m]
[m]
[m]
[m/year]
[mm/year]
[m/year] [mm/year]
[m/year]
[mm/year]
7 GANO
450554.732
5431524.174
701.17
0.00049
0.35
–0.00067
0.28
0.0014
2.16
9 HAGA
427582.955
5454992.094
1530.04
–0.00094
0.45
5.00E–05
0.34
–0.00069
2.59
11 KACI
448794.152
5470266.647
553.68
–0.00102
0.38
0.00026
0.3
–0.0018
2.26
18 KAWI
425866.824
5453835.261
1983.40
–0.00074
0.36
0.00025
0.27
–0.00058
2.13
12 LIES
404054.365
5469346.882
692.45
–0.00046
0.35
0.00089
0.26
0.00017
2.04
13 LOMS
442665.23
5449447.723
2634.03 –8.00E–05
0.36
0.00055
0.27
–0.00184
2.16
6 ROHA
398341.867
5434026.968
736.68
–0.00076
0.33
0.00021
0.25
–0.00071
1.90
15 SDOM
438492.561
5445087.551
1675.03
–0.00077
0.46
0.0006
0.35
0.0009
2.69
16 SKPL
443989.232
5448586.154
1772.47
–0.00084
0.35
0.00059
0.27
0.00044
2.18
22 TAKO
450764.513
5451388.624
820.29
–0.00091
0.62
0.00037
0.45
0.00123
3.51
19 ZISE
408976.34
5449663.668
1916.18
–0.00121
0.41
0.0001
0.31
–0.00124
2.25