Abstract
Damavand Mountain is a large volcanic system located in the north of Iran. Damavand has not had any recorded activity for 1000 years, but evidences of fumarolic activity have been observed near the summit. Due to the lack of data, there are no geodetic studies of volcanoes in Iran. Damavand magma source is modelled in this paper using Global Positioning System (GPS) instruments and Bernese software was used to process and support the interpretations. The maximum height variation obtained from GPS data was 4.7 mm. Direct modelling carried out in this paper demonstrated that a sill-like source model is the most suitable model to reproduce Damavand surface deformation.
1. Introduction
Volcanoes and eruptions are dramatic surface manifestation of dynamic processes within the Earth. They are mostly localized along the boundaries of the tectonic plates. Anyone who has witnessed volcanic activity has approved the variety and complexity of visible eruptive phenomena.
Analysis of the data obtained from adequate earthquake monitoring system indicates that nearly all eruptions have been preceded and accompanied by measurable changes in the physical and chemical states of the volcanic system. While geochemical methodologies of volcano monitoring have proved increasingly sophisticated, seismic and geodetic (ground deformation) techniques remain the most widely used tools in volcanic surveillance.
Systematic measurement of ground deformation was initiated at a few active volcanoes in Japan and the USA in the early 20th century (Dzurisin Citation2006).
Ground deformation due to volcanic magma intrusion is recognized as an important precursor of eruptive activity at a volcano. Before an eruption occurs, the ground surface generally expands due to the pressure increase within the shallow magma chambers located several kilometres below the volcano (figure ).
Figure 1 Ground deformation caused by magmatic activity (Abidin Citation1998)
The pattern and rate of surface displacement and the depth and the rate of the pressure increase within the subterranean magma chamber give important information about the state of a volcano. As ground deformation tends to precede eruptions by periods of hours to month, geodetic monitoring is an effective tool for hazard mitigation.
Ground deformation monitoring techniques for volcanic environment have developed from precise spirit levelling to electronic distance measurement (EDM) techniques and more recently to the use of Interferometric synthetic aperture radar (InSAR) and campaign-style or continuous GPS surveys (Janssen Citation2007).
Using elastic theory and data from several case studies, Mogi (Citation1958) developed a mathematical model in order to determine the position and depth of magma source based on ground deformation measurements and their locations. He assumed that the ground deformation is caused by a spherical source located below a volcano edifice which exerts hydrostatic pressure upward to the ground surface (Mogi Citation1958).
Volcanic deformation sources include inflating, deflating and growing bodies of various shapes and sizes which are collectively known as volumetric sources. Observed surface deformation can be fit to the prediction of source models. The modelling of source deformation, however, does not provide a unique description of the source causing the deformation (Dzurisin Citation2006).
Damavand stratovolcano is a volcano system located on the faulted/folded basement of central Alborz at approximately 5670 m above mean sea level. Damavand is about 10,000 years old and is classified as an active volcano (Omidian & Oskooi Citation2008). Shirzaei et al. have presented the first InSAR deformation, time series at the dormant Damavand volcano in northern Iran, over the period A.D. 2003–2008. The high-resolution data show a lateral extension of the volcano at the relative rate of ∼6 mm/year accompanied by subsidence at the rate of ∼5 mm/year at the volcano summit. They found that lateral motion of the east flank is more significant than that of the west flank (Shirzaei et al. Citation2011). Modelling of Damavand magma source using GPS data obtained in 2006 and 2007 is discussed in this paper.
2. Models of magma's deformation sources
Three mathematical source models are described in the following sections.
2.1. The elastic half-space: an immediate approximation of the Earth
In the mathematical source models discussed herein, the Earth crust is modelled as an ideal semi-infinite elastic body known as an elastic half-space. The half-space has one planar surface bounding a continuum which extends infinitely in all directions. The half-space is assumed to be made up of homogeneous and isotropic material. Hook's law specifies a linear relation between displacement at any point in the body and applied forces are, therefore, applicable to the elastic half-space source model (Dzurisin Citation2006).
2.2. Point pressure source
The point pressure or point dilatation source is often called the Mogi model. Kioo Magi has concluded that geodetically measured elevation changes and horizontal displacements associated with eruptions in Japan and Hawaii are resulted from inflation and deflation of magma bodies within the volcanoes (Mogi Citation1958).
The surface displacements with hydrostatic pressure changes within a limited spherical cavity embedded in an elastic half-space. The radius of the half-space is much smaller than its depth(α << d), hence:
is the radial distance from the centre of the cavity to a point on the free space. Equation (Equation2.1) is often further simplified by assuming υ = 0.25. The scaling coefficient includes the pressure change (ΔP) in cavity, its radius, shear modules (G) and Poisson's ratio (ν) of the half-space.
The displacements are approximately symmetric with the vertical displacement peaking directly above the source, while horizontal displacement obtains its maximum value at 
(figure ) (Mindlin Citation1936; Dzurisin Citation2006).
Figure 2 Profiles of axisymmetric displacement (vertical (red), horizontal (blue), magnitude of total displacement (black)) (Dzurisin Citation2006)
The surface displacement vector is radial to source (figure ) and its magnitude varies with the inverse square of the distance from the centre of the buried cavity:
Figure 3 Surface-displacement vectors and their components from inflation of a point pressure source, located below the origin at a depth d. All displacements are axisymmetric and the surface displacement (black arrow) is radial to the centre of the buried source (Dzurisin Citation2006)
Substituting the measured displacements would conclude the depth of the estimated source. Also assuming some values for source radius (α) and elastic constants (υ) and (G) the pressure increment (ΔP) would produce the deformation.
After a pressure increment (ΔP) on the inner surface of a spherical cavity, the cavity radius would be increased by (Δα), where
2.3. Sill-like magma chamber
Sill-like magma intrusions or chambers are represented by finite rectangular tensile dislocations, pressurized oblate spheroids or finite pressurized horizontal circular cracks.
A horizontal disk tensile dislocation is assumed in this study to illustrate the characteristics of surface deformation produced by a deep pressurized sill-like magma body embedded in an elastic half-space (figure ).
Figure 4 Pressurized degenerate oblate spheroid (horizontal point crack) used to represent a sill that is deep relative to its radius (Dzurisin Citation2006)
The disk tensile dislocation corresponds to a degenerate oblate spheroid. The length of vertical axis is much less than the length of horizontal axis. Although this model is approximate, it is adequate when the radius of the sill is much less than its depth. The free surface displacements are given by:
Figure 5 Profiles of horizontal (blue) and vertical (red) displacement compared with those from an equivalent spherical pressure source. The grey profile which is near to the red profile shows the Mogi vertical displacement and another grey profile shows Mogi horizontal displacement. The displacements are normalized by the power of the equivalent spherical source multiplied by the inverse square of the depth (α3ΔP/4Gd 2) and the distance is in source depths (Dzurisin Citation2006)
The zone of deformation is smaller, with the maximum horizontal displacement occurring at the distance of 
from the centre of chamber (Okada Citation1992; Dzurisin Citation2006).
3. Data and analysis
GPS data from the Iranian Permanent GPS Network for Geodynamics (IPGN) are used in this study. The IPGN covers all over Iran including Damavand. Speed and strain fields used for the stations were obtained in 2006 and 2007. Magma sources models introduced in the Section 2 of this document are correlated to surface deformation; thereby shape of Damavand magma source is hypothesized.
Six stations the positions of which are shown in figure have been used. Due to lack of data related to the summit of the mountain, the data related to stations around the summit for day of 200 in 2006 and 2007 were used in this study.
Figure 6 Position of stations around the summit
Bernese software (version 5.0) was used to process the data. Since baselines length was between 10 and 100 km, they have been considered as the baselines with middle length and quasi-ionosphere free (QIF) method was used for phase ambiguity resolution. In this method L1 and L2 should be processed similarly.
Baselines were produced by two methods which are star method and maximum observation. In the star method, baselines are produced through connecting a reference station to others. Reference station might be selected manually or automatically, so that in both cases, sum of baselines length are minimized. This minimization is an important factor for phase ambiguity resolution.
Figure illustrates baselines selected by the star method. Results from GPS data processing using star baselines in a geocentric reference frame are shown in .
Table 1 Coordinates resulting from star method for 2006 (upper table) and 2007 (lower table)
Figure 7 Baselines resulting from star method
In the maximum observation method, baselines are considered according to the number of common observation for their stations.
From all possible compositions, a set of baselines with maximum common observations is selected. This strategy guarantees that maximum number of observations from network could be used in double-differences processing. Figure shows the baselines resulting from the maximum observation method.
Figure 8 Baselines resulting from maximum observation method
The results of processing with this method are shown in . A comparison between standard deviations resulting from the two methods is given in and . According to proximity of final coordinates, the GPS data processing is concluded to have been carried out correctly.
Table 2 Coordinates are resulting from maximum observation method 2006 (upper table) and 2007 (lower table)
Velocity vectors for these stations were calculated during 2006–2007. Calculated coordinates, velocity (horizontal) and precision for each station as input for GMT software are used. Velocity vector with error ellipse is shown in figure . The amount of elevation changes (mm) is also indicated in figure . The resulting velocities with associated error ellipses are presented in and .
Table 3 Velocity components related to Damavand stations 2006–2007
Table 4 Semi-major axis and relevant azimuth of error ellipses (2-sigma errors) of Damavand stations
Figure 9 Velocity vectors and elevation changes related to stations 2006–2007
The stations illustrated in figure indicate a movement from southwest to northeast.
The length of baselines introduced between stations shows an average shortening of about 4 mm which is consistent with the results from geodynamic network of Iran from 2005 to 2007. In fact, Jamur el al. (2008) demonstrated that there is 4.3 mm shortening in baseline lengths in eastern Alborz (Damavand) per year.
Because of the exertion of different forces (like magma movement from a source toward the surface) on a deformable body, some changes on its shape and position occur, which happens slowly or suddenly. In order to monitor Damavand's deformation, in this case, strain fields were calculated in two dimensions. Using Delaunay finite element method, the stations were triangulated and strain tensors for the centre of gravity of each triangle were computed. In the Delaunay finite element method, triangular elements between a set of points are configured, so that each triangle would be close to an equilateral one. In the formation of each element, the assumption of homogeneous material is still retained. This assumption approves the harmony between magnitudes and directions of all obtained strains (figure ).
Figure 10 Strain fields around Damavand
Figure shows the strain field resulting for Damavand. Calculated strain field indicates a contraction that is compatible with the obtained results from the regional geodynamic network of Iran. Both cases display amount of 10−8 strain per year. Average standard deviations for strain components are in the order of 10−9 strain per year.
It is obvious that in large explosive eruption, a great volume of lava comes up to the surface of the earth causing changes in volume and pressure in volcano magma chamber. Due to lack of data related to hydrostatic pressure or volume changes in Damavand magma source, it is not possible to suggest a model with reliable information. Studies on other volcanoes around the world show that Campi Flegrei volcano in Italy is similar to Damavand in terms of geological characteristics so we applied the same models for Damavand.
3.1. Point pressure source of Damavand
Point pressure (Mogi) source is a current model for many of volcanoes for inverting surface deformation information of volcanoes surface (elevation changes) to source parameters (hydrostatic pressure change, depth and radius).
Normal surface deformation (uzs ) forms because of a volume change (ΔV) in Mogi.
Source which is given by Gottsmann et al. (Citation2006):
Assuming 
, a radius equal to 0.7 km, volume change equal to 0.5 cubic kilometre (as for Campi Flegrei) and various depths for magma chamber, elevation changes for 2006–2007, are presented in figure .
Figure 11 Elevation changes resulted from GPS and Mogi source with different depths
As it is inferred from figure , elevation changes from a source with 2 km depth and 0.5 cubic kilometre volume change are consistent with elevation changes resulting from GPS data.
Figure represents a diagram prepared for a chamber with various volume changes and 2 km depth. Figure indicates that the elevation changes from a source with 0.5 cubic kilometre volume change and 2 km depth fits better with GPS data.
Figure 12 Elevation changes for stations from GPS data and a Mogi source with different volume changes
Finally, figure shows the elevation changes from different depths and volume changes.
Figure 13 Elevation changes resulted from GPS data and a Mogi source with different volume changes and depths
As figure shows, a Mogi-type source with 2 km depth and 0.5 cubic metre volume changes fits better with GPS data. It should be noted that, generally resulting elevation changes do not show a uniform agreement with elevation changes resulted from GPS. In the case of four stations (ABSD, POOL, TEHN and BLDH), final elevation changes are close but it is not the same for two other stations.
3.2. Sill-like magma chamber of Damavand
There are some elastic solutions to explain deformations related to sill-like magma source presented. It is expected that surrounding cliffs behave elastically when magma comes up.
In this study, it is assumed that the deformation is due to an elastic half-space and that the source is under a uniform pressure. Maximum elevation change imputable a sill-like magma source is calculated by (Fialko et al. Citation2001):
Symbols mean the same as previous equations. Assuming 
, a depth equal to 3 km, pressure change equal to 4.5 MPa (as for Campi Flegrei) and changing radius for magma chamber, resulting elevation changes for 2006–2007 are shown in figure .
Figure 14 Maximum elevation changes for stations using GPS and a sill-like source with different radius
Figure shows that a magma chamber with a radius of 0.5 km is in more agreement with the GPS results.
To consider a source with 0.5 km radius and a pressure change of 4.5 MPa, figure is drawn for different depths.
Figure 15 Maximum elevation changes for stations using GPS and a sill-like source with different depths
Figure demonstrates that a magma chamber with 3 km depth fits better with GPS data.
Figure shows the elevation changes resulted from a source with 3 km depth, 0.5 km radius and different pressure changes. Because this model assumes that the radius of the chamber is less than the depth, depth of 3 km is considered.
Figure 16 Maximum elevation changes for stations using GPS and a sill-like source with different pressure changes
Figure shows that a magma chamber with a pressure change of 3 MPa fits better with GPS data.
Generally, sill-like magma chamber with 0.5 km radius, 3 km depth and 3 MPa pressure change matches better with the GPS data, even when compared to the Mogi source model. The maximum elevation change obtained from sill-like chamber is 4.6 mm compared to the maximum elevation change of GPS at 4.7 mm.
4. Conclusions and recommendations
Can Damavand become active? According to the evidence, it was activated 7300 years ago, but other evidences point to activity between 2000 and 3000 years ago. Based on volcanological research, volcanoes less than 10,000 years old are considered active. There is a need for supplementary studies concentrated on the physical and chemical properties of spa springs, changes in depth and volume of magma sources, or geophysical studies like magnetic surveys that can detect magma movement, and particularly inflation of the earth surface.
The models presented herein indicate shortening about 4 mm in baselines between stations. The results illustrate very little difference between the four stations (POOL, BLDH, TEHN, ABSD). This little difference is a good reason to choose a sill-like model for Damavand volcano. However, significant differences between the two stations (PLOR, PLZI) require dedicated geological studies. The significant differences can be considered as a new subject in future research using geological methods. A sill-like chamber gives better results than a point source, but only more geodetic data can give us a better insight into the inflation and other processes of Damavand.
References
- Abidin , H Z . 1998 . “ Monitoring the volcano deformation using repeated GPS surveys: experiences and plans ” . In Proceedings of International Workshop on Advances in GPS deformation monitoring , 18 – 27 . Perth , , Australia : Curtin University .
- Dzurisin , D . 2006 . Volcano deformation , Heidelberg , Berlin : Springer . xiii-296
- Fialko , Y , Khzan , Y and Simons , M . 2001 . Deformation due to a pressurized horizontal circular crack in an elastic half-space, with applications to volcano geodesy . Geophys J Int , 146 : 181 – 190 . doi: 10.1046/j.1365-246X.2001.00452.x
- Gottsmann , J , Rymer , H and Berrino , G . 2006 . Unrest at the Campi Flegrei caldera (Italy): a critical evaluation of source parameters from geodetic data inversion . J Volcanol Geoth Res , 150 : 132 – 145 . doi: 10.1016/j.jvolgeores.2005.07.002
- Janssen , V . 2007 . Volcano monitoring deformation using GPS . J Spatial Sci , 52 ( 1 ) : 41 – 54 . doi: 10.1080/14498596.2007.9635099
- Johnson , D J , Sigmundsson , F and Delaney , P T . 2000 . Comment on ‘volume of magma accumulation or withdrawal estimated from surface uplift or subsidence, with application to the 1960 collapse of Kilauea volcano’ by P. T. DELANEY, D. F. McTigue . Bull Volcano , 61 : 491 – 493 . doi: 10.1007/s004450050006
- Mindlin , R D . 1936 . Force at a point in the interior of a semi-infinite solid . J Appl Phys , 8 : 195 – 202. .
- Mogi , K . 1958 . Relation between eruptions of various volcanoes and the deformation of the ground surface around them . Bull Earthq Res Inst , 36 : 99 – 134 .
- Okada , Y . 1992 . International deformation due to shear and tensile faults in a half-space . Bull Seismol Sos Am , 75 : 1018 – 1040 .
- Omidian , O and Oskooi , B . 2008 . “ Hazard estimation of probable eruption of Damavand volcano ” . In Proceedings of the 3rd Disaster Management Conference , Tehran , , Iran : University of Tehran . 2008 December 24
- Shirzaei , M , Walter , T R , Nankali , H R and Holohan , E P . 2011 . Gravity-driven deformation of Damavand volcano, Iran, detected through InSAR time series . Geology , 39 : 251 – 254. . doi: 10.1130/G31779.1