Estimating volcanic deformation source parameters with a finite element inversion: The 2001–2002 unrest at Cotopaxi volcano, Ecuador

Deformation at Cotopaxi was observed between 2001 and 2002 along with recorded seismicity beneath the northeast (NE) flank, despite the fact that the last eruption occurred in 1942. We use electronic distance meter deformation data along with the patterns of recorded seismicity to constrain the cause of this unrest episode. To solve for the optimum deformation source parameters we employ inverse finite element (FE) models that account for material heterogeneities and surface topography. For a range of source shapes the models converge on a shallow reservoir beneath the southwest (SW) flank. The individual best fit model is a small oblate‐shaped source, approximately 4–5 km beneath the summit, with a volume increase of roughly 20 × 106 m3. This SW source location contrasts with the NE seismicity locations. Subsequently, further FE models that additionally account for temperature‐dependent viscoelasticity are used to reconcile the deformation and seismicity simultaneously. Comparisons of elastic and viscous timescales allude to aseismic pressurization of a small magma reservoir in the SW. Seismicity in the NE is then explained through a mechanism of fluid migration from the SW to the NE along fault systems. We extend our analyses to further show that if future unrest crises are accompanied by measurable seismicity around the deformation source, this could indicate a higher magma supply rate and increased likelihood of a forthcoming eruption.


Introduction
Magmatic intrusions often cause deformation at the Earth's surface as a result of the corresponding stress and strain being transferred through the crust [e.g., Gudmundsson, 2012]. If the rocks through which the magma propagates are relatively cold and brittle, and the associated strain rates are sufficiently high, then the intrusion can be accompanied by measurable seismicity. On the other hand, if the rocks are preheated, weakened, and ductile, or the strain rates are comparably low, then the magma can migrate aseismically. By combining the measured deformation with the recorded seismicity patterns, or lack thereof, we can build up an understanding of how magma is supplied, stored, and transported beneath active volcanoes.
Cotopaxi (5897 m) is an active stratovolcano located in the Eastern Cordillera of the Ecuadorian Andes ( Figure 1). Despite being in a regional compressional setting, magma emplacement and the development of Cotopaxi has been enhanced by a horizontal least principal stress and vertical NNE-SSW striking faults in a transfer fault zone [Fiorini and Tibaldi, 2012]. Unlike other volcanoes within the Eastern Cordillera (e.g., Tungarahua, Sangay, and Cayambe), Cotopaxi displays a bimodal magmatic evolution with a rhyolitic and andesitic eruptive history and very limited evidence of intermingling [Hall and Mothes, 2008]. To the west is the more dacitic Western Cordillera, with the densely populated Inter-Andean Valley sitting between them . The symmetrical conal edifice of Cotopaxi is topped by a 1 km 3 glacier covering a 21 km 2 area [Jordan, 1983]. Given the increasing population density on the surrounding flanks and the presence of a large glacier, the lahar risk is high [Barberi et al., 1995;Mothes et al., 1998;Aguilera et al., 2004]. As many as 100,000 people reside in the path of a lahar, the same size as one that occurred in 1877 [Pistolesi et al., 2013].
A stratigraphic study of eruptive products from 1140 to the present found no obvious patterns in magma composition, eruptive magnitude, or periods of quiescence; renewed magmatic activity could thus take the form of a single explosive event or mark the beginning of a swarm of eruptions [Pistolesi et al., 2011]. Such clusters of events are thought to be initiated by an increase in magma supply rate from depth, perhaps modulated by a deeper reservoir that recharges a shallower system [Pistolesi et al., 2011]. Petrological and geochemical results also indicate that Cotopaxi's volcanic system undergoes frequent recharge via andesitic magmas [Garrison et al., 2011]. Whether this magma supply proceeds seismically or aseismically, and the implications that may have for the given supply rate and thus the risk of eruption, remains to be seen.
Since seismic monitoring of Cotopaxi began, long-period (LP), very-long-period (VLP), and volcano-tectonic (VT) events have been detected [Ruiz et al., 1998;Kumagai et al., 2010Kumagai et al., , 2013, as well as tremor associated with lahars [Kumagai et al., 2009]. LP events recorded between 1989 and 1997, in a 3 km wide column stretching from +4 to −8 km depth, coincided with the occurrence of VT events but were not related to any volcanic unrest. The LPs were thought to be caused by the interaction of a hydrothermal system (fed by glacier melt) and hot material at depth (without fresh magma supply), while the VTs were a result of the increased stress caused by water and steam in cracks [Ruiz et al., 1998].
An unrest period at Cotopaxi began in January 2001 with increased LPs, followed by swarms of VTs in November 2001, and VLP/LPs in June 2002 [Molina et al., 2008]. Concurrent contraction of an electronic distance meter (EDM) line on the NE flank was also recorded during the 2001-2002 unrest and interpreted as flank inflation [Molina et al., 2008]. Waveform inversion of a single VLP event and particle motion analyses alluded to a volumetric change of a subvertical crack 2-3 km beneath the NE flank, but the location is not particularly well constrained given the network size and geometry (two stations on the north and northeast flanks) [Molina et al., 2008]. The combined results led to a conclusion of an intruded magmatic dyke as the cause of the recorded unrest [Molina et al., 2008].
In this paper we use the complete EDM data set to test the hypothesis that an intrusion in the NE caused the observed deformation. We show that the spatial deformation field is far wider than just the NE flank and use inverse finite element (FE) models that account for material heterogeneities and topographical complexities to search for optimal deformation source parameters. The results are combined with the seismic observations and further forward FE models that incorporate temperature-dependent effects to provide first-order insights on the architecture of Cotopaxi's magmatic plumbing system.

Data
The EDM network surrounding Cotopaxi was first established by the Instituto Geofísico Escuela Politécnica Nacional (IGEPN) in 1987, with the latest additions added in 2004. However, weather and access conditions prevent measurements being made at the same temporal rate across all base stations. For the [2001][2002] 10.1002/2014JB011731 unrest period, seven baselines provide usable data, with an unfortunate lack of coverage on the southern flanks ( Figure 1). To extract the values for this period from the time series we took the difference between the means for the periods before (1997)(1998)(1999)(2000) and after (2002)(2003)(2004)(2005) the unrest. All seven baselines display a contraction, thus indicating a general inflation of the edifice that is not restricted to the NE flank. Magnitudes of the contraction are between 1 and 8 cm, and the associated errors (2 : 0.8-2.5 cm) were calculated by a bootstrapping technique (Figure 1).

Numerical Model Setups
We use COMSOL Multiphysics (v4.4) to build and solve the finite element (FE) models in this study. The structural mechanics module is used in conjunction with an optimization procedure to solve the inverse problem. Subsequent forward FE models couple the structural mechanics to a heat transfer module to explore the effects of temperature-dependent mechanical properties.

Model Geometry
In order to account for sources offset from the center of the volcano and the nonaxisymmetric spatial deformation pattern, we use a full 3-D model geometry. Steep relief is known to alter modeled deformation patterns when compared to flat surfaces [e.g., Wadge, 1998, 2000;Trasatti et al., 2003;Champenois et al., 2014], so we include the real topography of Cotopaxi in our models via a digital elevation model (DEM) using 90 m SRTM data. This also allows a more accurate representation of the EDM lines in the FE model, as in reality they are recording deformation across three dimensions. The deformation source is represented as either a sphere, or oblate or prolate spheroid of varying eccentricity. Our final model geometry is 50 km × 50 km × 30 km in extent, and meshed with >165,000 tetrahedral elements, with a higher mesh density around the edifice and source regions ( Figure 2).
The inclusion of faults in the models were beyond the scope of this paper given the inverse approach selected and the lack of accurate data on their locations. Future work that can map their location and extent and has additional computational power should aspire to include them, given they can impact a surface deformation pattern if they are proximal to the deformation source [Folch and Gottsmann, 2006]. We also neglected to include the summit glacier as a separate geometric entity with its own mechanical properties as synthetic tests proved that the impact of this on the deformation field was negligible.

Solid Mechanics
Mechanical heterogeneities in the subsurface significantly alter stress and strain distributions when compared to the simpler, unrealistic, homogeneous representations that are often used in analytical models [Beauducel et al., 2004;Masterlark, 2007;Kinvig et al., 2009;Geyer and Gottsmann, 2010;Hautmann et al., 2010;Hickey et al., 2013]. Therefore, to include the effects of a mechanically heterogeneous subsurface, we use the IGEPN 1-D seismic velocity model from the RENSIG network to delineate layers of different mechanical properties [Font et al., 2013]. V P and V S velocities are converted to a Poisson's Ratio, , density, , and dynamic Young's Modulus, E, using the equations presented in Brocher [2005] (for all model parameters see Table 1): As the volcanic edifice is not accounted for in the velocity model, we use a dynamic Young's Modulus of 10 GPa (roughly equal to half of the layer immediately below) to account for further microheterogeneities, fracturing and alteration that can cause significant weakening. All the resultant values are consistent with those found in volcanic regions [Gudmundsson, 2011] ( Table 2).
The model boundary and initial conditions are adapted from the benchmarked approach outlined in Hickey and  to be applicable in a 3-D model: the top surface is free, the bottom surface is fixed, and the lateral surfaces have a roller condition ( Figure 2). An infinite element domain borders the bottom and lateral edges of the model geometry to negate any boundary affects impacting the deformation results  A pressurized spheroidal cavity is used to represent the deformation source, and an infinite element domain prevents any boundary affects from arising. More details can be found in Hickey and Gottsmann [2014]. (b) The meshed model configuration. Tetrahedral elements are used to mesh the model geometry, with smaller elements and a higher mesh density around the deformation source and volcanic edifice. (c) Thermomechanical model setup. Basal heat flux, surface temperature, and thermal insulation boundary conditions set up a background geothermal gradient, which is perturbed by a defined source temperature. The temperature distribution is used to calculate the viscosity distribution via equation (8) which is then used in a viscoelastic deformation model. An example temperature distribution and subsequent viscosity contours are displayed. Shear modulus Table 2 GPa K Bulk modulus Table 2 GPa Objective function equation (4)  of the interior. The model forcing is supplied by a boundary load applied uniformly to the surface of the source cavity walls and can be converted to a complementary volume change using the three-dimensional expansion of the source surfaces (also benchmarked against the solutions of Delaney and McTigue [1994]).

Optimization
To solve the inverse problem with a FE approach we use the Optimization module in COMSOL. This searches for the optimum location (longitude, latitude, and depth) and overpressure of the source to fit the EDM deformation data within its error and automatically rebuilds the mesh after each iteration. It works to minimize an objective function by varying a set of given design variables (source parameters) within their defined range. We specify the objective function, J, as follows:  The results of an initial optimization are used to inform a sequence of subsequent inversions, which use increasingly tighter constraints on the design variables, so as to build up nested parameter constraint grids and ensure the most robust final solution. This procedure is repeated for each source shape and size. The modeled EDM baselines, M i , are calculated by subtracting the three-dimensional straight line distance between the two points after the simulation from the distance before the simulation.

Thermomechanics
The observation of LPs caused by the interaction of a hydrothermal system and hot material at depth [Ruiz et al., 1998] indicates that a substantial amount of heat remains in the crust beneath Cotopaxi. Further heating is likely to be caused by any fresh input of new magma. Where rocks are significantly heated, they are unlikely to respond in an elastic manner; viscoelastic effects are thus more applicable [Dragoni and Magnanensi, 1989;Ranalli, 1995;Jellinek and DePaolo, 2003]. Temperature-dependent viscoelastic effects are known to impact the way stress and strain is partitioned in volcanic settings [e.g., Del Negro et al., 2009;Gregg et al., 2012;Gottsmann and Odbert, 2014], so we incorporate this in a series of forward FE models, guided by the best results obtained from the optimization procedure outlined above.
To include the influence of temperature-dependent viscoelasticity we adapt the approach outlined in Hickey and  for application in a 3-D model domain. This couples a heat transfer module, to calculate a steady-state temperature distribution, with the solid mechanics, where the temperature value is used to define the viscosity of the viscoelastic component across the model domain. The viscosity, , is calculated via where A d is the Dorn Parameter (a material constant), H is the activation energy, R is the universal gas constant, and T is temperature. For the boundary conditions, the surface is set to 273 K, the walls of the source are set to 1173 K, the bottom is given an inward heat flux of 0.09 W/m 2 , and thermal insulation is used on the lateral surfaces which results in a general background geothermal gradient of 30 K/km ( Figure 2). Elevated temperatures in the vicinity of the source result in reduced viscosities. This setup also requires that the elastic parameters for the solid mechanics are specified in terms of their shear modulus, G = E∕[2(1 + )], and bulk modulus,  Figure 3. Source parameter convergence. The optimization solver varies the source parameters to reduce the objective function. We use the results of an initial parameter grid (a, b, and c) to inform the parameter constraints of the next inversion (d, e, and f ). This process of nested parameter constraint grids ensures the final solution (g, h, and i) is as robust as possible, with the smallest possible objective function (or misfit).

Parameter Convergence and Source Location
The sequences of nested parameter constraint grids ensured that the final solution was as robust as possible: the results of the first search were used as the initial conditions for the next inversion, and so on. Evolution of the optimal source parameters can be traced through each iteration of the linked inversions ( Figure 3). This approach is justified by the continual decrease in the objective function, J, which serves as an indicator of modeled misfit. The sequence was stopped when J showed no further decrease between the first and last iterations of an individual inversion (e.g., Figure 3i). Figure 3 also shows that the number of iterations required in each inversion to converge upon a final solution decreased as the sequence of linked inversions progressed, further demonstrating the robustness of the solution.
All three source shapes converge on a solution beneath the SW-SSW of the edifice at depths of 1-2 km above sea level, regardless of size or spheroidal eccentricity (Figure 4). Oblate sources generally provide the best fit to the data, with prolate sources the worst (Table 3). In the case of some of the poorer fits to the data, the solver underperformed and got "stuck" varying the source location in one dimension. This is what creates the linear features that are visible in Figure 4 for the spherical and prolate sources. Also interesting to note is a second grouped population of spherical sources, slightly to the south of the main congregation of solutions (further evident in Figure 5b). Although the misfits for this second group of source locations are slightly higher, the linear alignment may allude to the SSW-NNE trending faults playing a role in the deformation recorded.
The convergence of the deformation sources in the SW contrasts strongly with the source proposed by Molina et al. [2008] in the NE to explain the seismicity. This mismatch is striking. To further explore this discrepancy and to counter the fact that the optimization solver does not always fully exploit the parameter constraint bounds supplied (i.e., if the tested parameters of one iteration provide a poorer fit to that of the parameters tested previously, the solver will "walk" in a different direction), we performed a series of extra tests in which the location of the deformation source was restricted to the NE. The results were poor, with the resultant modeled EDM lines only fitting a maximum of two of the observed EDM recordings (J > 0.006 m). Thus the most likely source to fit the observed deformation is SW-SSW of the edifice, as shown in

Best Fit Model
Of the parameters tested, the best fit model of this study fits the most EDM observations (six of the seven, Figure 6) with the smallest objective function (J ≈ 0.001 m). The source is oblate in shape with a = 0.25 km and b = 1.50 km. It is located SSW of the edifice (X = −78.447 +0.014  (Figure 4). The overpressure, ΔP = 80 +10 −5 MPa, translates to a volume change of ΔV ≈ 20 × 10 6 m 3 . Uncertainties on the source parameters were constrained by rerunning the inversion with the recorded EDM values replaced by their upper and lower error bounds. Figure 5 shows that the resultant best fit model parameters fit well within the most frequent source parameter distributions from all of the inversions conducted, while Figure 6 shows how the modeled EDM baselines converged to the final best fit solution as the inversion progressed through the different iterations.

Source Overpressure and Failure Criteria
The best fit modeled ΔP seems high (80 MPa) when compared to the usual estimates of rock tensile strength (T 0 ) that are often employed as maximum estimates of magma reservoir pressurization [Gudmundsson, 2006]. These range from 0.5 to 9 MPa when measured in situ [Gudmundsson, 2012], or up to 31 MPa when calculated from laboratory tests where T 0 = 6.88K C [Zhang, 2002], and K C is the tensile fracture toughness with a value of 3.5 ± 1 MPa m 1∕2 at magmatic temperatures [Smith et al., 2009]. In situ fracture toughness is considerably higher, particularly for stratovolcanoes [Gudmundsson, 2009] but has yet to be directly related to in situ tensile strength estimates.
There are a number of reasons why ΔP could be high. First, the inversion models employed an elastic host rock rheology to reduce the computational requirements, which is known to overestimate pressure requirements [Trasatti et al., 2003;Newman et al., 2006;Del Negro et al., 2009;Hickey et al., 2013]. Linked to this are the values used to represent the elastic medium. We converted seismic velocities to a Young's Modulus and  Figure 5. Source parameter frequency distributions. The colored lines (blue, red, or green, for spherical-, prolate-, or oblate-shaped sources, respectively) represent the frequency histogram data for each source shape. All four parameters solved for in the inversions are displayed: longitude (a), latitude (b), depth (c ), and pressure (d). The source parameters used to define the best fit model are displayed with the diamond on the x axis. In Figure 5b, the "I" represents a second smaller congregation of spherical source locations albeit with a higher misfit than the main population of source locations.
a Poisson's Ratio and thus obtain dynamic values, but static values are perhaps more appropriate given the slower nature of volcanic deformation in comparison to the propagation of seismic waves [Gudmundsson, 2011]. As static moduli can be up to 13 times smaller than dynamic moduli [Gudmundsson, 2011], this would give a smaller ΔP estimate proportional to that factor for a linearly elastic domain. In this case, our modeled overpressure may actually be below any static reservoir failure criterion.
In reality, however, the rocks encasing a magmatic reservoir are unlikely to behave in an elastic manner due to high thermal gradients and numerous microfractures. Viscoelastic behavior is more likely in this case, which brings the applicability of a static limiting overpressure into disregard due to viscous dissipation of stress acting to relieve the building pressure and expansion of the reservoir boundaries by solid-state creep [e.g., Jellinek and DePaolo, 2003]. Hence, a dynamic failure criterion based on a critical strain rate is perhaps more appropriate [e.g., Gottsmann and Odbert, 2014], and the 80 MPa value for ΔP is of less significance.
Regardless of any line of argument presented above, there is also the possibility that in the present case the small reservoir SSW of the edifice did indeed "fail" and contribute to a migration of fluids to the NE. This is discussed in a later section.

Influence of Thermomechanics
When the effect of temperature is included in the models via a temperature-dependent viscoelastic setup, there is a reduction in the required overpressure to fit the same deformation pattern. The thermomechanical model setup was applied to the best fit model from the inversions (described above), and for the given settings ΔP was reduced by 25% to 60 MPa when applied as a constant throughout the simulation. The resultant ΔV remained at ∼20 × 10 6 m 3 (due to viscous expansion following the instantaneous elastic inflation) and explains how the same amplitudes of surface deformation can still be met.
The thermomechanical setup also introduces a time-dependent aspect into the model. However, as the temporal resolution of the EDM data is poor, and also varies between baselines, it does not warrant a full in-depth study of the influence of temperature-dependent viscoelasticity and time-dependent pressurization on deformation predictions [e.g., Gottsmann and Odbert, 2014]. Instead, we use the setup described to examine the rheology of the magma reservoir host rock. The settings employed in this study predict rocks with a viscosity of 10 14.3 Pa s immediately at the magma reservoir boundary (Figure 7). This result is dependent on the steady state temperature distribution that is calculated in the model (Figure 7), which of course may not reflect true subsurface conditions. A steady state temperature distribution requires a certain amount of time to reach equilibrium, so if this time has not passed, the heat from the hot zone(s) will not have conducted as far through the crust, and the general temperature distribution will be reduced. On the other hand, if a hydrothermal system is present (something the current model does not incorporate) heat could be convected away from the source more efficiently and the temperature distribution would be increased. Therefore, the temperature distributions presented here serve only as first-order estimates and assume a thermal legacy in the system resulting from the previous eruptive history.

Journal of Geophysical Research: Solid Earth
Consequently, we test how variability in the parameters used to calculate temperature and viscosity can influence the resultant rheology of the host rock. The values of thermal conductivity and heat capacity are somewhat obsolete in the steady state simulation as they only affect the time it takes for heat to conduct and materials to warm up, respectively. Geologically plausible variations in the geothermal gradient (20-60 K/km), input as variations in the basal heat flux, still result in viscosities of 10 14.3 Pa s around the chamber (although they do alter the background temperature which would influence the resultant surface deformation). Changes to the surface temperature also do not affect the viscosity of the host rock. The most important controls are the source temperature, and the activation energy (used in equation (8)). Magmatic temperatures at Cotopaxi have been estimated between 743 and 1673 K [Garrison et al., 2011], which when translated into a source temperature in the model gives viscosities between 10 12.7 and 10 17.4 Pa s at the host rock interface. Activation energies between 105 and 165 KJ/mol [Meissner and Tanner, 1992] give viscosities between 10 13.7 and 10 16.3 Pa s. The combined results give a potential range in viscosity of 10 12.7 -10 17.4 Pa s for the magma reservoir host rock. However, in reality, it is possible that the scenarios with the lower range of viscosities would become part of a "capture front" [Marsh, 1996], given that a viscosity of 10 16 is often used as an indicator of a magmatic "mush" zone. Either way, the host rocks would certainly respond in a viscoelastic manner.

Resolving the Deformation and the Seismicity
With the observed seismicty in the NE and the origin of the deformation in the SW-SSW, there is an initial discrepancy when trying to resolve the two with the same source. Here we explore a potential mechanism to reconcile the deformation and seismicity simultaneously.  (8). High temperatures around the source produce rocks that behave in an increasingly viscoelastic manner.

Aseismic Inflation
Given there was no elevated seismicity recorded in the SW during this unrest period, perhaps the processes that lead to the best fit reservoir being pressurized proceeded aseismically. This can be the case when rocks behave in a ductile manner, which is more likely when the rocks are presumed to be heated, weakened, and viscoelastic [Jaeger et al., 2007]. Ductility is then expected when the timescale of a forcing process is significantly longer than the Maxwell viscoelastic relaxation time of the rock [Dragoni, 1993].
To test the hypothesis that the inflation of the small magma reservoir inferred beneath Cotopaxi proceeded aseismically, we compare the elastic pressurization timescale to the viscoelastic timescale using the equations presented in Jellinek and DePaolo [2003]. The elastic timescale equation was derived for a spherical source but can still provide a first-order approximation when comparing the two timescales for our best fit oblate-shaped source. The elastic timescale, e , is where ΔP crit is the critical overpressure to fail a chamber, V ch is the volume of the chamber (in this case 4 ∕ 3 ab 2 ) and Q is the time-averaged magma influx [Jellinek and DePaolo, 2003]. The general Maxwell and aseismic deformation is expected if v << e . Hence, we can rewrite the equations in order to express them in terms of a required host rock viscosity value: To solve equation (11)  Magma is supplied to the reservoir inferred from the inversion of EDM geodetic data. The pressurization of this reservoir causes deformation at the surface and fluids to migrate (green arrows). The fluids preferentially travel along the NNE-SSW trending faults (dashed yellow lines) which produce the VLP and LP signals that are observed, while the VTs are a result of increased fluid pressure in cracks/pores in a colder, more brittle part of the crust. The best fit modeled deformation source (solid red ellipse) and topography are drawn to scale. See the text for more details. viscosity is 10 16.1 -10 18.6 Pa s. In the most part the broad range in host rock viscosity values inferred from the thermomechanical models (10 12.7 -10 17.4 Pa s) are indeed smaller, and we can thus propose that the reservoir pressurization occurred aseismically.

Fluid Migration and Seismicity
The seismicity source proposed by Molina et al. [2008] cannot explain the observed deformation. However, if we assume the NE seismicity locations to be robust, we can explain both the seismicity and the deformation with the best fit model inferred from the EDM measurements.
We propose a mechanism whereby magma supply and pressurization of the magma reservoir initiated a migration of hydrothermal and/or magmatic fluids. This could have been caused by one, or a combination of, a change in the local stress field, a temperature increase, or fresh fluid injection. The fluids would flow most efficiently along the NNE-SSW trending faults [Fiorini and Tibaldi, 2012], which could take them into the NE quadrant of the volcano from the source in the SW. Fluid movement and mass transport have long been proposed as the causes of LP and VLP seismicity [Chouet, 2003;McNutt, 2005] so this could explain the presence of both types of signals in the NE during the 2001-2002 unrest period. Where the fluids move away from a hot magmatic zone into colder parts of the crust, the rocks would become increasingly brittle. Hence, the VT signals could be the result of excess fluid and gas pressure in cracks and pores. Our mechanism, similar to that proposed for LPs and VTs at Cotopaxi between 1989 and 1997 [Ruiz et al., 1998], is summarized in Figure 8.

Implications for Eruption Forecasting
Clusters of eruptions are known to occur at Cotopaxi when there is an increase in the time-averaged magma supply rate, which may be modulated by a deep magmatic reservoir [Pistolesi et al., 2011]. We have inferred a small shallow reservoir from the inversion of EDM data, which could allude toward some connection between the two levels of magma storage. A column of hot material is proposed to exist beneath the volcano from depths of −8 km to +4 km which contributes toward the creation of LP signals [Ruiz et al., 1998]. This column may in fact represent a complex plumbing system that connects a deeper reservoir to shallower storage regions. Assuming this plumbing system geometry, we can infer that the current magma supply rate is low due to the lack of eruptions and aseismic inflation (following equation (9), a higher magma supply rate would decrease e and hence increase the likelihood of seismicity). If the magma supply rate increases and there is recorded seismicity accompanying the deformation, then this may be indicative of a critical level of unrest that precedes a forthcoming eruption. Work to constrain the deeper system through evaluation of a wider deformation field and the nature of the fluid movement through gravimetric studies should form part of a subsequent study.

Conclusions
Using the recorded patterns of seismicity and deformation, we have presented a viable mechanism to explain a period of unrest at Cotopaxi volcano between 2001 and 2002. A thorough examination of the EDM deformation data set indicated a general expansion of the entire edifice. Inverse finite element models that account for material heterogeneities and surface topography while searching for the optimum deformation parameters converge on a solution for a small, shallow source SW-SSW of the edifice. Sequences of nested parameter constraint grids ensured a robust solution. The best fit model is an oblate source with a = 0.25 km and b = 1.50 km, at the location X = −78.447 • , Y = −0.724 • , Z = 1.17 km, and a volume change of ΔV ≈ 20 × 10 6 m 3 .
A deformation source in the SW contrasts strongly with the source proposed in the NE by Molina et al. [2008] to explain the seismicity. However, predicted EDM baseline results were especially poor when the deformation source was restricted to the NE. Therefore, to reconcile the deformation and seismicity simultaneously, we use thermomechanical viscoelastic finite element models and comparisons of elastic and viscous timescales to show that the pressurization of the shallow magma reservoir proceeded aseismically. Subsequent fluid migration from the SW into the NE along NNE-SSW trending faults caused the seismicity due to mass transport and excess pore pressures. The lack of eruption following this unrest period and the aseismic nature of the magmatic recharge indicate that the magma supply rate was low. Future unrest episodes where measurable seismicity accompanies the deforming magma reservoir could be indicative of a higher magma supply rate and a critical level of unrest that could precede a forthcoming eruption.