Modeling 3D Free-geometry Volumetric Sources Associated to Geological and Anthropogenic Hazards from Space and Terrestrial Geodetic Data

: Recent decades have shown an explosion in the quantity and quality of geodetic data, mainly space-based geodetic data, that are being applied to geological and anthropogenic hazards. This has produced the need for new approaches for analyzing, modeling and interpreting these geodetic data. Typically, modeling of deformation and gravity changes follows an inverse approach using analytical or numerical solutions, where normally regular geometries (point sources, disks, prolate or oblate spheroids, etc.) are assumed at the initial stages and the inversion is carried out in a linear context. Here we review an original methodology for the simultaneous, nonlinear inversion of gravity changes and / or surface deformation (measured with di ﬀ erent techniques) to determine 3D (three-dimensional) bodies, without any a priori assumption about their geometries, embedded into an elastic or poroelastic medium. Such a fully nonlinear inversion has led to interesting results in volcanic environments and in the study of water tables variation due to its exploitation. This methodology can be used to invert geodetic remote sensing data or terrestrial data alone, or in combination.


Introduction
Over the past few decades, an explosion in the quantity and quality of geodetic data has been increasingly applied to geological and anthropogenic hazards [1][2][3][4][5][6][7]. It includes ground-, air-and space-based observations that span a range of spatial and temporal resolutions. The development of space geodetic techniques has been particularly important to the measurement of surface displacements, resulting in a sharp increase in the number of active areas where deformation has been detected [5][6][7][8][9]. These new sensitivities and capabilities have produced the need for new approaches for analyzing, modeling and interpreting geodetic data [4,5].
Simultaneously, the rapid growth of the global population concentrated in large cities has led to significant increase in risks when such highly populated areas are also subject to natural (e.g., flood, volcanic and seismic activity) and anthropogenic hazards [4][5][6][7][10][11][12][13][14]. Improving our understanding of the processes associated with these geological and anthropogenic hazards is crucial for developing systems that monitor and evaluate them, as well as to aid decision makers during a crisis.
In the case of volcanic hazards, detecting the accumulation and ascent of magma can provide advance warning of future eruptive activity [15]. Such magmatism is often associated with ground deformation and gravity change, making geodesy a basic tool for the monitoring of volcanic hazards [2,4]. Volcano geodetic observations also provide valuable data for exploring the geometry and volume of magma plumbing systems, which offer a critical context for interpreting the mechanisms and characteristics of unrest and eruption [16,17].
Additionally, land subsidence, ranging from local collapse to the broad regional lowering of Earth's surface, represents the main geomechanical effect related to the removal of subsurface support [5]. Subsidence can occur as a result of (i) natural factors (e.g., tectonic activity, self-consolidation of recent sedimentary deposits, oxidation and shrinkage of organic soils) [18,19], and (ii) anthropogenic processes (e.g., groundwater pumping [20][21][22][23], urban development [24], hydrocarbon or mining exploitation [25,26]). Land subsidence associated with overexploitation of aquifers is a hazard that commonly affects large areas worldwide.
In both cases, volcanic activity and subsidence associated with aquifer exploitation, the link between data and interpretation is given by direct models and inversion techniques, such that both are critical components to any geodetic monitoring study.
Typically, modeling of deformation and gravity change follows an inverse approach using analytical solutions for volume change or shear faulting within an elastic half space-computationally efficient methods that provide insight into the depth, geometry, and strength of subsurface magma bodies and slipping faults [27,28]. Normally, regular geometries (point sources, disks, prolate or oblate spheroids, etc.) are assumed at the initial stages [27] and the resulting inversion is carried out in a linear context.
In a volcanic context, surface displacements are inverted to infer valuable constraints on the active magmatic sources [2,17,29,30]. Surface deformation is a direct consequence of the dynamics of volcanic plumbing systems, and reflect the shape of magma intrusions, the volume of intruding/arising magma, and the emplacement mechanisms.
Advances in computing power have facilitated rapid development of numerical methods that are capable of both incorporating material and structural heterogeneity [31,32] and investigating the characteristics of the subsurface [33]. They can even be used in inverse approaches [34], blurring the lines between the traditional "inverse" and "numerical" categorizations [5].
Models that combine a suite of geodetic observations can provide exceptional insight into the properties of magmatic systems. Joint assessment of, for example, tilt, GNSS, and differential interferometric synthetic aperture radar (DInSAR) data, leverages the strengths of different data types (in this case, temporal resolution, 3D displacements, and spatial resolution, respectively) to better image magma storage and transport pathways [4].
Some of the most compelling combinations of various types of geodetic data involve joint modeling of deformation and gravity change [35,36]. Typically, gravity and vertical deformation are expected to correlate (accumulating magma results in both mass and volume increases, and vice versa), but a growing body of literature documents instances of gravity change in the absence of deformation [36,37]. Such results may be caused by magma filling or draining from void space [38], magma compressibility [39,40], source geometry and depth [40], or topographic effects [41]. In addition, providing a window into the characteristics of magma plumbing systems, joint modeling of deformation and gravity change can constrain the physical properties of the magma. Similar (and sometime the same) models can be used [5] to study other kind of geological and anthropogenic hazards, such as land subsidence related to groundwater pumping. This represents a hazard commonly affecting large areas worldwide, often associated with the increasing demand upon groundwater resources due to expanding metropolitan and agricultural areas in semiarid and arid regions [21]. The surface ground deformation thus constitutes a signature of the processes in the reservoir and can provide information about those subsurface processes [42]. Many recent studies have focused on this topic.
The modelling of surface deformation patterns can provide significant insights into the temporal changes of pore pressure as well as the 3D geometry of a reservoir in response to its exploitation over the time [43,44]. A number of different techniques have been developed in recent decades to estimate the surface deformation pattern related to volume changes in elastic and poroelastic media [23,[44][45][46][47][48][49][50]. Inverse modeling is required to achieve success in such an endeavor [1,41]. A proper understanding of the subsidence mechanism is essential to calibrate protocols and best practices for monitoring natural and anthropogenic phenomena, with the aim to reduce vulnerability and risk for infrastructures, economies, natural environments and human life.
Analytical models are used to estimate the amplitude and pattern of surface deformation based on assumptions about the media and perturbation source (e.g., using elastic or poroelastic theory) [17,50,51]. They provide a relatively simple method to model surface deformation for reservoirs of any geometric shape. Furthermore, given that these techniques assume that most of the surface deformation is explained by the poroelastic expansion or contraction of the reservoir, less in situ geological data is required than that needed for numerical models [50].
Usually, analytical magma source models represent mathematical abstractions based on strong assumptions (e.g., a priori definition of the source shape, use of geometrically simple geometries, and selection of medium characteristics) that make the set of differential equations describing the problem tractable. Sometimes the analytical approach must be limited to a very simplified and rough approximation of the real source. Nevertheless, analytical models may be useful in certain applications where priority is given to the fast computation of inverse solutions.
On the other hand, numerical Finite Element (FE) modeling allows realistic features such as topography and crustal heterogeneities to be included. However, FE modeling requires large computational resources and it is still very time-expensive to solve the inverse problem with FE models for real time monitoring purposes, i.e., it can take several hours of computation to invert data in order to estimate the source parameters [31,52].
In an operative real time monitoring system, minimizing the delay in the estimation of location and characteristics of the shallow magma source and their evolution in time is the desired goal. Minimum delay in inferring the possible location and timing of an impending eruption would provide critical time for decision makers. Indeed, the fast response of a magma source model can be crucial in early warning for mitigating the possible dramatic consequences of an eruption near populated zones. In this regard, the use of an analytical model for data inversion could speed up the source estimation. Nevertheless, it is still necessary to impose assumptions on the a priori source geometry. Unfortunately, this selection often is highly unconstrained due to the complexity and lack of knowledge of volcano plumbing systems.
Camacho et al. [1] developed an original methodology aimed at the determination of the 3D geometry and location of the causative bodies by inverting ground deformations and gravity changes due to pressure and/or mass anomalies embedded into an elastic medium. Such a fully nonlinear inversion has led to interesting results in volcanic environments, where ground deformations are related to over-pressured magmatic bodies [17,53,54] or water extraction [5]. This methodology can be used to invert geodetic remote sensing data or terrestrial data alone, or in combination.
In this work, we present a review of this method for the simultaneous, nonlinear inversion of gravity changes and surface deformation using bodies with a free geometry. Assuming simple homogenous elastic or poroelastic conditions, the method determines a general geometrical configuration of pressure and density sources. These conditions do not completely represent the complex geological conditions of active areas but are do obtain good results as described in the literature. These sources are described as an aggregate of pressure and density point sources, fitting the entire data set (given some regularity conditions). The approach works in a step-by-step growth process that allows us to build very general geometrical configurations.

System of Nonlinear Equations
Camacho et al. [1] and Fernández et al. [5] consider a set of observation points at the surface (e.g., a geodetic network composed of several stations) (x l , y l , z l ), l = 1, . . . , n, where gravity changes dg i , i = 1, . . . , n g , and position changes (dx j , dy j , dz j ), j = 1, . . . , n p , have been observed, nearly simultaneously, but not necessarily in the same locations.
Let us denote by d the n g +n p observation vector of components dg i , dx j , dy j , dz j . They try to model this vector by using m buried point sources, located in (X k , Y k , Z k ), k = 1, . . . , m, and characterized by corresponding values for volume v k , (positive or negative) incremental pressure p k , and (positive or negative) incremental density ρ k changes.
They fit the data d as where d c represents the n g +n p vector of modeled values dg c i , dx c j , dy c j , and dz c j for gravity and position changes, and ε is the n g +n p vector for residual values coming from inaccuracies in the observation process and also from insufficient model fit [1].
They simultaneously calculate the modeled changes, dg c i , dx c j , dy c j , and dz c j , by means of simple formulas for vertical attraction and deformation effects that are due to additional mass and pressure of several point sources (X k , Y k , Z k ; v k , ρ k , p k ) within a homogeneous elastic halfspace characterized by a shear modulus µ (given in stress units of pascals, Pa) and a Poisson's ratio of ν 0.25.
The surface deformation dx, dy, dz, due to a buried point source is considered mainly as composed by (1) the effects that are due to the incremental pressure p k and expansion radius within the elastic half-space, described by the Mogi model [31], and (2) the effects produced by the loading of the additional point mass v k × ρ k within the elastic half-space, described by the Boussinesq model [55,56].
The surface gravity changes dg can be decomposed as follows [1]: 1.
Free-air effects, corresponding to the relocation of the benchmarks, due to elevation changes according a free-air vertical gravity gradient (about −290 µGal/m for Campi Flegrei); this effect can be included in the model fit using modeled or observed elevation changes; 2.
Newtonian effects due to density changes within the original boundaries of the deep bodies; 3.
Newtonian effects due to mass relocation or change of volume; 4.
Effects due to mass uplift in the surface corresponding to elevation changes. These effects can be obtained using another vertical gravity gradient, depending on the regional terrain density (similar to the Bouguer correction [57][58][59][60]). Effects 1 and 4 depend on the surface elevation changes and can be combined by using a combined gravity gradient F (about −210 µGal/m); 5.
Water table effects which correspond to very local and shallow perturbations; 6.
Topographical effects used for the approximate approach of assuming different source depths at each point [61,62].
Furthermore, they assume possible common and constant terms dx 0 , dy 0 , dz 0 , for the position changes (for instance, an unknown change of the leveling origin) and a possible constant term, dg 0 , to represent the simplest form of a long-wave component, a global offset value, or a common uncorrected groundwater effect. Then, with the conditions described previously, they adopt the following expressions for the modeled changes [1]: where G is the gravitation constant, ν, µ are elastic parameters (Poisson ratio and shear modulus), and g 0 is a mean surface gravity value (g 0 ≈ 9.8 Gal). Equations (1)-(5) represent linear equations for ρ k and p k , but a nonlinear relationship with respect to the geometrical parameters (point sources "filled" with non-zero anomalous density or/and non-zero anomalous pressure). We observe that they calculate the free-air and "Bouguer" effects of the elevation changes as proportional to the modeled values dz c j , instead of proportional to the observed ones. This allows us to include gravity and elevation changes corresponding to different stations and gives better final results.
In the case of DInSAR data, deformation ds on the surface points is measured along the radar line of sight (LOS) corresponding to the direction of the satellite antenna. Then, the modeled changes ds c i for the n data points can be written as where α, β are the direction angles (azimuth and incidence) for the antenna pointing direction.
The described system is valid for a volcanic loading problem [1], but if we want to model displacements produced by water extraction in an aquifer [5] we have to move from an elastic to a poroelastic medium. Then, considering the linear theory of poroelasticity [63], the horizontal and vertical components (du, dv, dw) of the movements at a point (X, Y, Z) of the free surface, due to a differential nucleus, located at (x, y, z) and with sides ∆x, ∆y, ∆z, corresponding to the reservoir with local overpressure ∆p, are [1,45,46]: where ν denotes Poisson's ratio (normally estimated at~0.25) and c m is the uniaxial compaction coefficient. Assuming that displacements at the surface are almost directly proportional to the thickness ∆z of the reservoir, the volume integrations for a parallelepiped cell of sides ∆x, ∆y, ∆z and overpressure ∆p in equation (2) can be simplified to integration in the horizontal plane only given rise to [5,46]: where: The integrals I i for the displacements in the i-direction are: I x (p, q, r) = arcsinh p I y (p, q, r) = arcsinh q p 2 + r 2 (12) This formulation provides the direct calculation of the surface effect of a single parallelepiped cell. The total effect of an anomalous structure described as aggregation of m small parallelepiped cells is obtained, according [45], as the summation of the partial effects. This direct formulation can be used to carry out the inverse approach ito determine the 3D pressure source structure responsible of the observed surface deformations for the aquifer case.

Misfit Conditions
Camacho et al. [1] assume a Gaussian distribution of the observation uncertainty, described by a covariance matrix Q D for the gravity and position data dg i and dx j , dy j , dz j (or ds j ), then the minimization condition for residuals ε, ε T Q −1 D ε = min, leads to the maximum likelihood solution. It is also assumed that the data are not correlated, and Q D is a diagonal matrix of estimated variances corresponding to observation of gravity and position changes. Gravity and elevation observations are independent and, then, Q D can be decomposed into two covariance matrices Q G (n g × n g ) for gravity and Q P (n p × n p ) for position changes. Then the minimization condition for observation residuals can be written as where ε G and ε P are n g and n p vectors for gravity and elevation residuals and θ is a weighting factor for the balance between gravity and deformation fits. This factor is introduced to allow a more versatile fit, making it possible to give optional priority to one kind of observable data over the other. Values for θ can be adopted according to the accuracies of the data sets and considering the approximate equivalence of 1 cm ≈3 µGal.
As is typical for the inversion of geophysical data, problems of singularity and instability for the solution can arise because of inadequate data coverage (e.g., often the number of data points is smaller than the number of unknowns), because of inaccuracy of the data, and because of intrinsic ambiguity of the problem (i.e., if we assume that positive and negative anomalous density and pressure can be contemporaneously present in the model). A possible approach for solving these problems is to use fuzzy logic [64] based on a general estimation of the properties of the physical environment. A more classical and simple process to avoid the instabilities of the solutions is obtained by means of additional minimization or smoothing conditions for the norm of the solution model as where the model vector m is constituted by the values ρ k , p k , k = 1, . . . , m, for the cells of the model and Q M is a suitable covariance matrix corresponding to the physical configuration of cells and benchmarks. This matrix provides a balanced model, avoiding very shallow solutions. We use a normalizing diagonal matrix Q M with elements q k , k = 1, . . . , m, given for volumes ν k and distances d jk as Condition (14) is a stabilizing term [65,66] to control the incremental mass and the pressure for the bodies (weighting according to matrix Q M ). It prevents very large fictitious values of mass and/or pressure resulting from a rather poorly determined model (for instance, one that is due to the coupling of some positive and negative sources, aligned stations, peripheral sources).
A mixed minimization equation, is finally adopted as the constraining Equation (1) for residuals and for model magnitude. λis a factor for selected balance between fitness and smoothness of the model. Low λvalues produce a very good data fit (even with noisy inversion) but cam produce an irregular and/or extended model. Conversely, high λ values produce concentrated and smooth models but often with a poor data fit, potentially with autocorrelated residuals. The optimal choice is determined by an autocorrelation analysis of the residual values, as that value producing a null (planar) autocorrelation distribution. See [67] for more details.

Exploration Approach for Solving the System
The model system, (1)-(5) or (8)- (12), must be satisfied within the minimization constraint (16). The system constitutes a nonlinear problem of optimization with respect to the geometrical properties. Iterative or explorative methods are mostly used for this purpose [68]. Iterative methods (gradient methods) are supported by a long tradition, but they generally require starting with a good initial solution. The more versatile explorative methods can be advantageously applied when the model space to be explored is reasonably sized (small number and narrow range of parameters). In this case, random explorative processes (e.g., genetic algorithm or simulated annealing) are frequently adopted. However, taking into account the very large number of degrees of freedom required to describe the density and pressure models as aggregations of thousands of small cells filled with anomalous values, a general exploratory inversion approach applied simultaneously for all the cells, would be ineffective. An alternative approach is to build the anomalous 3-D structures by means of a growth process and apply an exploratory approach for each step of the growth process.
We carry out a step-by-step process of growth of the 3-D models, using an exploratory technique to find a new cell to be filled with density and/or pressure and aggregated to the models. Therefore, for the kth step of the growth process, k cells have been filled with the prescribed anomalous values for pressure and density, giving rise to modeled values d c . For the new (k + 1)th step, we try to find and fill a new cell to fit the system: where f > 1 is a scale factor to allow for a fit between the anomaly of the provisional (not totally developed) model and the observed anomaly (gravity and position changes). To solve, we calculate the value e 2 = ε T Q −1 D ε + λ f 2 m T Q −1 M m for each of the empty cells according to an exploratory technique. In the (k + 1)th step we choose, as an optimal cell to be filled, the jth cell, one giving The process continues until we reach a scale factor value f 1. In the final result, we arrive, nearly automatically, at 3-D models described as the aggregation of point sources filled with the prescribed anomalous values for density and pressure. These simultaneously fit the observed values for gravity and position changes and maintain a small size for both anomalous models.
Camacho et al. [55,69] and Gottsmann et al. [70] provide some simulation examples showing the suitability of this 3-D inversion approach as applied to gravity data. Results are quite satisfactory. A trend towards producing rounded bodies or other distortions is observed for zones of low-reliability, largely because of the additional condition necessary to solve the underdetermined system. Camacho et al. [1] and Fernández et al. [5] carried out some simulation examples to study the performance of the previously described inversion process, divided into two groups. The first group represents volcanic sources as a sum of pressure and mass sources and the second represents water table variation as pressure change sources. It is observed that the fit between the simulated and modeled structures is quite satisfactory [1,5]. The corresponding inversion structure fits the location, depth, and magnitude, but offers a more rounded geometry because of the imposed smoothing condition.
In addition, the sizes and depth are quite representative, although some distortions are observed for very peripheral or very deep areas in the model.

Application Results
The interpretation methodology, as described above, has been applied to different and complementary test cases. Some of these are volcanic, focusing on the study of a particular period, a long-term time series, or the real time inversion of deformation data acquired during a volcanic crisis. An additional application is the study of the deformation associated with water extraction to determine the volume of the extracted water.
These test cases show the applicability of our methodology for the interpretation of geodetic remote sensing data alone or combined with other terrestrial ones, using only displacement or deformation data, or in combination with gravity variations.
The obtained results for the interpretation methodology demonstrate their applicability to the study of a variety of natural and anthropogenic hazards with different approaches. We describe the application test cases in the following subsections.

Modeling of Campi Flegrei Unrest 1992-2000 Using Deformation and Gravity Changes
Camacho et al. [1] applied the methodology to actual data from the caldera of Campi Flegrei (Italy), Figure 1. They used gravity and leveling data for terrestrial gravity benchmarks and LOS deformation observed by DInSAR for ascending and descending passes. All the data correspond to nearly the same time period, 1992-2000. They employed change rates, i.e., mean annual displacement or gravity changes. See [1] for a description of the geological setting, the volcanic history, and previous geophysical models. The leveling benchmarks employed are only a part of a longer network covering the whole area [71][72][73], referenced to an Istituto Geografico Militare Italiano (IGMI) station located in Naples. The double-run survey meets first-order standards. Uncertainties in elevation differences are less than ±2 mm D 1/2 , where D is the distance along the line between benchmarks [70], and a final uncertainty of the elevation value d z at each benchmark result is typically less than 1 cm [74].
In the time interval detained in this work, the gravity network consisted of 15 stations (Figure 1a), each one coinciding with benchmarks of the leveling network. Gravity-height data have an inverse linear correlation [73,75]. As an example, Figure 1b shows the temporal changes for gravity and elevation at the Serapeo station, since it is one of the oldest stations and is close to the area of maximum vertical movements. Camacho et al. [1] observe that the elevation change data show a generally continuous feature in time, but the gravity data contain a higher noise level, with a larger oscillating feature. Then, instead of the real observed data (the values just observed for 2000 and 1992, independent of the other values), they use some filtered (or interpolated or extrapolated) values. This approach allows them to avoid local or instantaneous perturbations or errors in the data, providing a better definition of the inversion model. In addition, they use a mean deformation rate (mGal or centimeters per year) for each station as input data for the inverse approach. See [1] for additional details on leveling and gravity.
For the same area, Campi Flegrei, Small Baseline Subset (SBAS) results (linear velocity and displacement time series) were obtained from a set of 165 European Remote Sensing Satellites 1 and 2 (ERS-1 and ERS-2), and 62 ENVISAT SAR data acquired between 1992 and 2008 on ascending (track 129, frame 809) and descending (track 36, frame 2781) orbits [76]. Camacho et al. [1] use a temporal subset of this data set overlapped with the common time span of the leveling and microgravity data. In particular, data were limited to the periods 1992-2000 for the descending pass and 1993-2000 for the ascending pass.
Ascending and descending deformation linear rate maps are shown in Figure 2. Both maps show a roughly circular pattern of subsidence (increase in line-of-sight phase change) with maximum average velocities larger than 3 cm/yr. Despite some minor uplift [1], the linear velocity during the whole period 1992-2000 is negative. Linear ground deformation rate and time series evolution have been extensively validated against independent geodetic techniques [77]. The global precisions of the SBAS results were determined to be ±0.5 cm and ±0.1 cm/yr for the linear ground deformation rate and time series displacement evolution, respectively [78]. See Camacho et al. [1] for more details in the DInSAR data sets. In summary, the input data sets used [1] are: (1) The terrestrial data set in terms of linear deformation velocity, including 15 benchmarks with coordinates (x, y, z), gravity change rates dg, and elevation change rates dz for the period 1992-2000; and (2) the total DInSAR data (approximately 15,000 points). For a faster process they selected pixels separated by a distance larger than, for instance, 300 m. It produces two subsets, each with approximately 1330 data points. They apply the inversion methodology to the combined subsampled data set (about 2670 data points) with both ascending and descending data and the terrestrial gravity and elevation, using all the data simultaneously.
A suitable relative weighting for the different types of data is necessary to carry out a simultaneous fit [1]. Leveling data have a good confidence level (better than 1 cm), and their temporal evolution is quite stable (e.g., Figure 1). Gravity data also have a good confidence level (better than 10 mGal), but their temporal evolution is noisier (Figure 1). For the DInSAR data, we have a rather good accuracy level (assumed to be ±0.5 cm and ±0.1 cm/yr). In accordance with these considerations, we assume a similar relative weighting factor for our various types of data (gravity, leveling and DInSAR). As a result, the modeling process is going to be conditioned by these weighting criteria.
The assumed values for the elastic parameters are Poisson ratio = 0.25 and shear modulus = 10 GPa [79]. The assumed gravity gradient used is −220 mGal/m. It corresponds to a free-air value of −290 mGal/m [72] and a value 70 mGal/m for the Bouguer correction for Campi Flegrei, resulting from the analyses carried out by Berrino et al. [80].
To obtain the inverse models they consider a 3-D partition of the local subsurface volume into 14,500 elemental point sources (located at the center of the cells with sides of about 120 m). Using this distribution of "empty" elements, the inversion approach produced the models for pressure and density as defined by the aggregation of filled point sources according to the prescribed extreme values for anomalous density and pressure. After some trials, they selected the extreme values of ±20 kg/m 3 and ±10 MPa as a suitable contrast to obtain extended anomalous bodies [1]. The standard deviation of the inversion residuals is 0.2 cm/yr for the elevation changes, 0.6 mGal/yr for the gravity changes, and about 1.1 cm/yr for the LOS values. Figure 3 shows the 3D source model obtained, and Figure 4 shows some additional views of the inversion residuals for the case of the ascending LOS component.

Discussion
The first main result is that the process does not produce any significant volume for anomalous positive or negative density. The resulting anomalous model for density changes consists of only a very few filled cells located close to the surface below some benchmarks, which correspond to very local effects. The anomalous 3-D depressurized model ( Figure 3) [1] shows a clear structure located below Pozzuoli at a mean depth of about 1500 m bsl. The vertical profiles in Figure 3 suggest a "partially filled parabolic glass" shape along a WNW-ESE course, with a less-developed similar shape along the orthogonal course. The picture also suggests some small tracks from the main body toward the surface.
In this application case [1], the anomalous mass model obtained from a joint inversion of gravity and leveling changes is nearly empty. There is no anomalous mass below the central area. Only some very shallow (depth about 200 m) and small anomalous (positive and negative) masses are adjusted close to some benchmarks, corresponding to local effects. It suggests that mass changes are not detected below the main deformation area. The gravity changes observed in the benchmarks of the main deformation area can be fully interpreted as produced by free air and Bouguer effects of the uplift (according the assumed gravity gradient). They [1] conclude that the phenomena below the central area do not involve input or output of significant mass (magma or fluids).
The absence of a significant mass change would imply no temporal variability of the fluids in the source region [1], which might be consistent with an exponential decay of fluid pore pressures during the analyzed slowdown subsidence phase (post-1987), following the rapid uplift phases (1982)(1983)(1984) during which high pore-pressure conditions were established. A physical mechanism that could reconcile these observations could be a long-lasting transient pulse of pore-pressure diffusion confined to the caldera filled material, which does not necessarily imply a mass transfer process (fluids in and out) in the hydrothermal system [81]. However, the absence of mass change does not preclude high pore-pressure conditions from producing fluid flow at the more permeable, shallow, hydrothermal system levels, as suggested by the geochemical variation of volatile discharge at Solfatara [82]. Indeed, local fluid flow could explain the few very shallow mass change sources required to fit some benchmarks.
All deformation data (terrestrial, ascending, and descending DInSAR) provide a similar low-pressure structure located at a shallow depth (about 1500 m) below the survey area [1]. The shape of that structure nearly fills the bottom of a parabolic cup with the same additional filling in the walls ( Figure 3). Therefore, they did not detect signs of the deep magmatic source for this period. Some main WNW-ESE elongation is detected at about N100 • E. The parabolic cup is superimposed with depth on the shape of some existing structure, namely the inner caldera edge, as depicted by a 2.5-D inversion of on-land and offshore gravity data [83]. Moreover, its bottom corresponds to the depth of the bottom of the inner caldera (between 2 and 2.5 km). The modeled structure fits well the polygenic body inside the caldera. We interpret this as corresponding to the dynamics of the shallow (depth 1-2 km) hydrothermal system, as previously suggested by several authors [74,84,85]. The hydrothermal system is confined to the filling caldera materials and is limited by the inner caldera structure.

Modeling of Campi Flegrei Unrest 1993-2013 Using Only Displacement Data
DInSAR is used extensively today for mapping ground deformation with high spatial resolution and sub-centimeter precision over large areas and is a suitable tool for ground deformation monitoring of active volcanic areas [4,53]. Spatial resolution of modern synthetic aperture radar (SAR) sensors ranges from 1 to 20 m over areas from 10 × 10 km 2 to 250 × 250 km 2 [53]. For modern satellite constellations the repeat cycle ranges from 1 day to a few weeks, but the typical repeat cycle for a single satellite mission ranges from 12 to 41 days. Repeatedly acquired SAR data from a single sensor are used to obtain line-of-sight (LOS) time series of ground deformation by applying small baseline subset (SBAS) [86][87][88], persistent scatterer [89] methods, or a combination [90]. The results of these techniques, however, are limited to the time period of the particular data set and do not distinguish horizontal and vertical motion.
The multidimensional SBAS (MSBAS) technique [53,91] combines multiple DInSAR data sets into a single solution with improved characteristics: Lower noise, denser temporal resolution than earlier DInSAR time series, and continuous temporal coverage. The MSBAS methodology is a natural extension of the original SBAS method that addresses the data redundancy and multidimensionality of the problem by decomposing LOS DInSAR measurements into vertical and horizontal (mainly east-west) time series of surface deformation using ascending and descending DInSAR data. MSBAS has been applied to the mapping of anthropogenic [92,93] and natural [91] ground deformation, successfully producing 2-D time series with dense temporal resolution and high precision.
Samsonov et al. [53] apply the MSBAS DInSAR processing technique to 20 years of ERS-1/2, Envisat, and RADARSAT-2 data from Campi Flegrei. The resulting spatiotemporal signals were determined with a duration and resolution that had not been achieved before in this or any other volcanic region. The used radar data are listed in Table 1. The theoretical derivation of the MSBAS technique is described in detail in [53,91]. MSBAS processing starts after both ascending and descending images are acquired [53], so the time series begin on the first date of the ascending ERS data (10 January 1993). The resulting 1271 highly coherent interferograms were used in the MSBAS processing, producing vertical and east-west time series spanning 20 years (1993-2013) and 385 time steps. The average error in the east-west time series is approximately 0.07 cm, while the error in the vertical time series is approximately 0.09 cm, although this varies spatially and can be as much as twice that in regions of maximum deformation (supporting information in [53]). In regions such as Campi Flegrei, where there is a significant trend, this estimate provides the worst case scenario, and the actual error is likely smaller. See [53] for more details on the DInSAR processing. Figure 5 presents the results of the MSBAS processing [53]. Figure 5a shows the vertical change in surface height between the initial and final time steps. The associated deformation time series is shown in Figure 5c for one location and presents a more complicated picture than the net subsidence of Figure 5a,b (deformation time series for selected locations are presented in the supporting information by [53]).  (Figure 5c). Additional times series can be seen in [53].
Using this new data set, Samsonov et al. [53] determine the potential sources below the caldera using the described methodology [1]. They apply the inversion method to specific time periods related to the inflation and deflation at Campi Flegrei observed in Figure 5c, 1993-2013. Table 2 shows those time periods and the resulting parameters for the associated inflation or deflation sources. The various periods of subsidence and uplift (deflation and inflation) can be modeled as separate sources. The results are summarized in Table 2 and in the next Figures. We show the results from the nonlinear inversion method for different time periods (1993-1999, 1999-2000, 2000-2005, 2005-2007 and 2007-2013) and the actual and modeled surface deformation rates (cm/yr) for both vertical and east-west motion for each of the considered time intervals together with the residuals between the observed and modeled displacements. See Figures 6-15. The models resolve the actual displacements quite well, with the exception of coherent residuals in the Solfatara region that are likely associated with the shallow source under the fumarole field [94] and not modeled in that work.  Images were saturated at the corresponding scales in order to improve comparison and highlight differences among panels [53].  Images were saturated at the corresponding scales in order to improve comparison and highlight differences among panels [53].  Images were saturated at the corresponding scales in order to improve comparison and highlight differences among panels [53].  Images were saturated at the corresponding scales in order to improve comparison and highlight differences among panels [53].  Images were saturated at the corresponding scales in order to improve comparison and highlight differences among panels [53].

Discussion
The computed deformation by Samsonov et al. [53] show consistent caldera subsidence for the period 1993-1999. However, deformation changed from subsidence in 2001 to primarily uplift in 2005, in conjunction with changes in geochemistry indicators such as CO2/H2O ratios [95], indicating an increase in relative enrichment of a magmatic CO2 component in sampled fumarolic gases. The rate of uplift continued to increase through the present, reaching more than 5 cm/yr.
Comparison of the deflation and inflation sources for the various activity periods shows that they are remarkably similar in location, shape, and volume (Table 2 and Figures 6-15). The inflation source tends to be deeper, nearer the bottom of the marine precaldera and postcaldera deposits and limited in depth by the transition to the gas-bearing deposits layer [53]. The deflation source is generally shallower, near the top of the marine precaldera and postcaldera deposits. This behavior suggests that inflation-deflation cycles are generated by sources fed from below and deeper than the 3 km detected here. Migration of fluids seems to be focused at the top of the gas-bearing formation, overcoming lithostatic pressures and effectively transferring stress to the surrounding rock during inflation periods. This is followed by gradual relaxation of the overpressure (e.g., pore fluid pressure diffusion) at shallower depths [1]. While others have modeled deformation at Campi Flegrei with two simple sources at 1700 m and 3600 m in depth [94], here the bottom of the inflation source extends below 3000 m. In this case, one extended source model with unconstrained geometry explains the deformation [53]. Samsonov et al. [53] showed that the deformation driving mechanism is strongly controlled by the caldera structure and stratigraphy and must involve large extended sources in a layered hydrothermal system. Due to the location of Campi Flegrei and its proximity to the city of Naples, more active monitoring for signs of possible imminent eruption is warranted.
Apart from the geological and geophysical conclusions, their study demonstrates that a temporal resolution on the order of one or two weeks can be achieved for MSBAS time series, nearing that of continuous GPS (cGPS) daily time series, by simultaneously processing DInSAR data from different sensors. The accuracy of MSBAS vertical and east-west components is very high, potentially higher than cGPS. At the same time, the spatial resolution of MSBAS results is superior to cGPS. Therefore, the MSBAS approach can be used for monitoring volcanic activities in near real time, particularly in remote areas where other observation techniques are unavailable or insufficient.

Real Time Tracking of Magmatic Intrusions During Volcanic Crisis: 2008 May Etna Eruption
The goal followed by Cannavò et al. [17] is to estimate in real time, for a given volcano showing signs of unrest, the location, size and volume distribution of the deformation source and to track those parameters through time without any a priori assumption on the source geometry. They seek to delimit the area with the highest likelihood of being the location of an opening fissure or vent for the impending eruption. They apply the approach by Camacho et al. [1] to model volcanic sources and track their evolution in time. Their objective is to address the limitations in the inversion process dictated by the running time and the a priori assumption of the shape of the intrusion or magma chamber in order to produce a useful real time that only adds seconds to the time needed to obtain displacement data. This methodology [17] performs a nonlinear inversion of the displacement data, producing extended bodies [1,53] with a free geometry within seconds (less than a minute) in a fully automatic way, i.e., without requiring any user input.
To demonstrate the applicability and results obtained using this methodology in real time tracking of magmatic intrusions using ground deformation data they [17]  The intrusion produced large ground displacements captured by the volcano cGPS network [96].
Due to the magnitude of the event and the quality, quantity and variety of available data, this eruption was a good benchmark to test new volcano modeling techniques. In particular, this study models the ground deformation pattern associated to the co-intrusion dyke propagation in the shallower part of the northern sector of Mt. Etna. The high rate GPS (HRGPS) data for 10 stations (Figure 16) from the Etnean GPS permanent network (Etn@net) is employed to illuminate the dynamics of the magma intrusion [97]. A real time scenario is simulated by computing HRGPS solutions from the raw data and running the inversion code, thus simulating a real processing and inversion stream during the eruption process. The real time analysis starts the day before the eruption when GPS signals showed some transients, though it should be noted that the shallow co-intrusion period lasted from the 8.30 am UTC [96] (when the seismic swarm started) on 13 May 2008 to 04.00 pm UTC the same day (when almost all the seismic activity had ended). The real effusion [93] phase, together with lava fountaining, started at 9.30 am UTC, as revealed by Meteosat Second Generation satellite. The fountaining stopped at around 11.00 am UTC [98]. This implies that the early warning window to localize the impending eruption site was only 1 hour (from 8.30 am to 9.30 am UTC).
Cannavò et al. [17] reproduce the real time scenario of the eruption. Though the test data are processed after data collection (post-processing), real time operation was simulated. Displacements are calculated by a robust linear regression [99] of the GPS solutions in the 30-minutes window. An interval of 30 minutes was chosen as a tradeoff between the minimum response time of the warning system and the minimum level of noise in the estimated displacements ( Figure 17). See [17] for more details on the GNSS processing method.  [100]. In green, the cumulated sources during the yearlong recharging phase and in orange the sources active during the day of the eruption. Cannavò et al. [17] suggested fast magma ascent from the reservoir level (2-3 km depth bsl) along paths indicated by the source geometries and the earthquake locations (dashed lines in the figures) [17].
Inversions are performed in a continuous way from the simulated real time displacement data computed every 15 minutes. There is some general uncertainty for any inversion approach with respect to the pressure contrast ∆P of the causative structure. Cannavò et al. [17] use 1 MPa as the value of pressure change for inversion in this test case, considering typical values obtained studying different eruptions for Mt. Etna volcano and from other studies of the considered test case eruption [96]. Authors perform several tests, considering different values of pressure change (0.5 to 3 MPa), that produced very similar results.

Results
Before testing the real time response of the inversion procedure, Cannavó et al. [17] use the inversion methodology to analyze the inflation phase preceding the 13 May 2008 eruption. They model the time evolution of the magmatic source for the period 1 January 2007 to 12 May 2008 by considering 15 day-intervals (Figures 16 and 18). It is determined that, during the months of June and July 2007, the system was subject to recharging episodes characterized by an ascending high-pressure body. As shown in Figure 19, the modeled body ascended from a depth of about 10 km, to reach a level of about 2-3 km depth (b.s.l.). This modeling study allowed estimation of the rate of magma ascent and suggests an approximate 3D geometry. It was assumed that an a priori estimated knowledge and interpretation of mid-term (weeks-to-months) ground deformation signals are routinely obtained in most volcano observatories. Such estimated models account for typical monitoring resources, which include a weekly report of volcanic activity, integrating the available geodetic data on those time scales.
Nine months after the deep magma ascent, a dyke intruded in the North-East sector of the Etna summit area. The simulated real time inversion produced a time series of 3D body sources shown in Figure 20. In particular, they observed that in the evening (20-22 UTC, all the times are in UTC hereinafter) of the day before the eruption, a vertically elongated source from 2 km bsl to 2 km asl (above sea level) identifies a magma batch likely migrating towards the surface (see the Supplementary Videos [17]). Starting from 6.00 am of the day of eruption, a large positive pressure source located at about 2.5 km asl in the NE sector of the summit area triggered the dyke intrusion at 7.00 am towards the shallow complex source in the NW-SE direction ( Figure 19, panel a). The latter source reached its maximum expansion in that direction at 11.30 am (Figure 19, panel c) and contracted only after 12:30 pm. This activity exactly follows the sequence of events and models reported in previous studies [96,98]. The estimated pressure sources show, minute by minute, that most of the recorded deformation corresponds to a shallow dilation process located close to the NE sector of the summit area, just where the eruption fracture opened.  Figure 16 shows the traces of the major sources that emerged from the inversion for both the studied time scales (long and short) under consideration together with seismicity during these periods. The uncertainty of the source location is estimated at about 500 m in the horizontal component and 900 m in the vertical component [17]. Thus, by considering the associated errors of the earthquakes (200-300 meters in all the 3D components) recorded by the local permanent seismic network managed by INGV, from Figure 16 they can estimate the earthquake swarm to be co-located with the estimated magma sources during the intrusion phase. It is worth noting that all the identified elongated model sources seem to occur near the periphery of the High Vp Body (HVB) [100,101]. A comparison between the shape and the position of the volumes estimated by Cannavò et al. [17] and the ones modelled by Bonaccorso et al. [98] for the period 1993-2000 suggests that there is a magma reservoir at about 2 km bsl and a vertically elongated source located in the south-western border of the HVB, where magma migrates towards the surface. Figure 20 provides a 3D-view sketch summarizing the results.
The goodness-of-fit for the estimated time varying model is assessed by computing the root mean square error (RMSE) between the measured and the modelled displacements. Discrepancies between observed and modelled values are in the order of 0.2 and 0.4 mm/day for horizontal and vertical components displacements, respectively [17].  Some spurious sources appear scattered around in the domain (see e.g., Figure 19 and the Supplementary Videos by [17]). They primarily are y due to the residuals in the data (non-ideal network topology configuration, noise measurements, and errors in the parameter settings) that the algorithm attempts to account for. A suitable trade-off between good data fit (with possible aggregation of spurious sources) and good model stability (with poor data fit and too simple model geometry) is obtained by analyzing the autocorrelation of the resulting data residuals [1]. Nevertheless, spurious sources are usually easily distinguishable, but sometimes they aggregate, forming a volume that may appear as a real source. In such cases, considering the spatial location, the size, the shape and their temporal evolution, they can usually be identified and discarded as true sources of deformation.
In order to quantify the reliability of the solutions obtained with this method, Cannavò et al. [17] carry out an uncertainty analysis. Figure 21a shows the decrease of the sensitivity to pressure with the depth and with the distance to the stations. In fact, the required pressure (for a 1 km 3 body) to produce a particular magnitude of surface deformation at the current network stations (root mean square, rms, 1 mm) increases with depth and distance to the stations. For deep and distant pressure elements in the inverse model, the relative uncertainty increases with a pattern similar to that of Figure 21a. Figure 21b,c show the 3D uncertainty distribution for vertical and horizontal displacements of the isolated pressure sources. It also is valid for some geometrical components or details of extended sources.
13 Figure 21. 3D uncertainty maps for the considered GPS network of 10 stations, representing, respectively, (a) the pressure changes, (b), the depth changes, and (c) the horizontal deviations required to produce a surface deformation with quadratic mean value 1 mm. Map (d) represents the minimum horizontal deviation for an isolated body of 1 km 3 with a pressure change given in panel (a) to produce a deformation at the GPS stations with rms magnitude of 1 mm [17].
These maps show the location values (in km), for a source volume 1 km 3 and a pressure value 1 MPa, required to produce a particular surface deformation level (1 mm, rms value). They can be interpreted as showing the pattern of relative uncertainty for the geometrical aspects (depth and horizontal location) for the inverse model. In addition, the displacement uncertainty pattern is heavily dependent on the pattern of the uncertainty for the source intensity (pressure and volume) of the assumed source element.
A more valuable uncertainty map for model geometrical aspects is obtained by assuming equally sensitive source elements (according to the pattern of Figure 21a). Figure 21d shows the sensitive 3D pattern for displacements for a volume of 1 km 3 , but assuming pressure values increasing with distance according to Figure 21a. These last model uncertainties show a different pattern, essentially due to geometrical aspects. Further details on uncertainty analysis can be found in the Supplementary material by [17].
The uncertainty analysis shows the power resolution of the GPS network for the proposed source estimations. The uncertainty maps show that the sources are well-resolved in the volumes where they have been estimated in Figures 16, 18 and 19, thus making the results more reliable.
To analyze the uncertainty of the results to the GPS network configuration, Cannavò et al. [17] carried out the same sensitivity tests reported in the manuscript by using the latest GPS network on Mt. Etna. Here they estimated the results for the current configuration of the GPS network composed of 40 stations. The results show a rather different pattern, mostly corresponding to the geometrical distribution of stations. For some common areas the 40 stations sensitivity is obviously higher (smaller pressure changes or smaller displacements) than that corresponding to 10 stations. This means that the more stations involved in the inflation/deflation process (i.e., measure displacements), the better the resolution of the model [17].

Discussion
The main goal of this test case is to study the ability of the described inversion approach to instantaneously provide a distribution of pressure sources for real time volcano monitoring purposes. Cannavò et al. [17] show that this inversion methodology overcomes the limitations of fixed-shape analytical models and the time-consuming nature of numerical approaches based on finite elements, providing a tool to infer the timing and location of the potential impending eruption, thus helping to activate warning procedures in good time. It may therefore be a crucial tool for civil protection purposes. For example, during the first hours of the studied 2008 Etna eruption, bad weather conditions prevented researchers from identifying the location of the ongoing eruption. However, the use of a tool like the one proposed here could have assisted researchers during the peak of the crisis by helping define the magma location in real time and infer the area of the opening eruptive fracture.
The hypothesis of homogeneous elastic medium does not appear to affect in any critical way the quality of the test results considered here, as was expected with the time periods used [103,104]. They show good agreement with previous results described in literature, e.g., the shallow NW-SE oriented source described above overlaps the previously modeled dyke [96].
Tracking the evolution of the magmatic source beneath the volcano surface in real time has always been a challenge for researchers working in modern volcanic surveillance. The proposed methodology is able to estimate, in real time, the location, size and shape of the active magma source in volcanic areas and can track those parameters through time, without any a priori assumption on the source geometry. This information can be the basis for forecasting activity in a monitoring and hazard assessment system to assist decision makers during a volcanic crisis. To this end, this test case showed the effectiveness of the proposed methodology for one of the best monitored active volcanoes on Earth. Mt. Etna represents a significant hazard to the local population living in neighboring areas and to buildings and infrastructure. The methodology, though successfully applied here to a Mt. Etna eruption, can be tuned and used in all active volcanic areas. Compared to other approaches in volcano monitoring by means of geodetic data, the advantages of the proposed tool are: (1) The speed of data inversion (few seconds); (2) the absence of assumptions on the shape of the source, (3) obtaining a 3D geometry of the magmatic sources and their time evolution; and (4) the easy implementation of the tool (it runs on a common personal computer without any special requirements).

Land Subsidence Associated with Overexploitation of Aquifers: Lorca, Spain
The Lorca region, located in the Alto Guadalentín Basin of southeastern Spain (Figure 22), is affected by subsidence rates of up to 10 cm/yr as a direct consequence of long-term aquifer exploitation [5,[21][22][23][24][25][26][27][28][29][30][31] ( Figure 22). This region is characterized by semi-arid climate conditions, with average precipitation rates of 150 mm/yr and an average annual temperature [105] of~18 • C. The basin is infilled with Quaternary alluvial fan systems overlapping Tertiary sediments transported by the Guadalentín River along the depression located in the eastern part of the Betic Mountain Range (an ENE-WSW oriented alpine orogenic belt resulting from the Nubia-Iberia ongoing convergence [106][107][108]).
The Guadalentín Basin aquifer is composed of two contiguous sub-basins: Alto and Bajo Guadalentín (Figure 22). From a hydrogeological point of view, the basement beneath the aquifer is composed of several relatively impermeable Paleozoic metamorphic complexes overlain by permeable Miocene conglomerate and/or calcarenite series. The top of the succession comprises Pliocene-Quaternary, low-permeability, compressible conglomerates, sand, silt, and clays [21,109]. The Alto Guadalentín aquifer covers an area of approximately 277 km 2 . Historically, piezometric levels were located closely to the land surface, allowing the development of a number of artesian wells and permanent lagoons [109].
Available piezometric information consists of water-level time series for a few points, from the 1970s to present. A long drought period from 1990 to 1995 (also in 1999-2000 and 2005-2007) reduced natural recharge and increased pumping in the Guadalentín Basin, which led to an increased resources deficit. All this information indicates a long-term trend in the consumption of groundwater resources [5,21].
Interferometric synthetic aperture radar (DInSAR) studies, while detecting the high subsidence rates affecting the Alto Guadalentín Basin, also identified a delayed transient of nonlinear compaction of the Alto Guadalentín aquifer due to the 1990-1995 drought period [21]. This suggested a relationship between local crustal unloading and stress change on active faults bordering the basin [51]. Later work [22] extends those studies using advanced differential DInSAR (A-DInSAR) techniques to process ALOS PALSAR (2007-2010) and COSMO-SkyMed (2011-2012) radar images. The combination of multi-sensor SAR images with different resolutions allowed for a longer monitoring time span of 20 years (1992-2012) over the Alto Guadalentín Basin. Additionally, the satellite measurements study area. Furthermore, the work presented a new soft soil thickness map and collected historical piezometric data, in order to assess aquifer system compressibility and groundwater level changes in the past 50 years. From the analysis of these data with A-DInSAR displacement measurements, the authors conclude that the governing mechanism of the Alto Guadalentín aquifer system is an inelastic, unrecoverable and delayed compaction process between water level depletion and ground surface displacement, related to the presence of very thick (>100 m) unconsolidated sediments (clay and silts).
Despite the aforementioned achievements, the previous studies focusing on the deformation in the area are based on DInSAR analysis using ascending and/or descending acquisitions, without any combination of the datasets to estimate both vertical and horizontal (E-W) components [21,22]. Therefore, only the line-of-sight (LOS) displacement field is known in the Alto Guadalentín area at a regional level and it was assumed to correspond completely to vertical displacement. Although this is a common procedure in subsidence studies using DInSAR measurements [112][113][114][115][116][117][118], the main consequences are (i) the neglecting of possible horizontal displacement components and (ii) the likely overestimation of vertical displacement.
Here, with the decomposition of LOS measurements in the E-W and vertical components over the investigated area, additional constraints are provided for the spatial and temporal evolution of the subsidence process as well as on the main governing mechanisms (e.g., temporal changes of pore pressure, geometry of the reservoir). With this primary aim, Fernández et al. [5] established a GNSS network consisting of 33 stations in 2015, which densely covers the Alto Guadalentín basin (Figure 22). This network was observed in survey mode. They analyze the measurements carried out in November 2015, June-July 2016 and February 2017 [5]. GNSS raw data was processed by adopting standard processing strategies for this type of network and referred to a local reference frame in order to estimate the 3D deformation field (see [5] for additional information). Despite the limited time interval covered by the surveys, they estimated, for the first time, a significant 3D deformation field, which is primarily related to the local exploitation of the aquifer. SAR data from the Sentinel-1 Copernicus constellation, acquired in ascending and descending orbits for the same time period, also were processed to obtain the respective LOS displacements. Using the GNSS and SAR-based deformation fields, both the vertical and horizontal components of the displacement over the entire area were estimated [5]. Finally, they use these results for comparison between them and interpretation using the forward model and inversion technique previously described.
Fernández et al. [5] carried out three Global Navigation Satellite System (GNSS) campaigns, in November 2015, June-July 2016 and February 2017. The 3D velocity field results are shown in Figure 23 for both the vertical and horizontal components determined by comparing the coordinates obtained for the time spans of the three surveys. See [5] for additional details on the observation and processing methodologies, and on the results. The maximum vertical subsidence rate (9.0 ± 0.5 cm/yr) is of the same order of magnitude as that previously detected by earlier DInSAR studies [21,22]. The maximum horizontal displacement rate detected is 2.5 ± 0.3 cm/yr (about 28% of the vertical displacement rate), a non-negligible amplitude. In the area showing the highest subsidence rate, again previously detected by DInSAR techniques, a characteristic pattern of horizontal deformation appears ( Figure 23). These deformations, as theoretically expected, show a centripetal pattern towards the zone of maximum subsidence, located in the central part of the monitored area.

15
In the southern area, where there is a relative maximum in the LOS displacement detected by DInSAR, significant horizontal motions are detected ( Figure 23) associated with GNSS stations 23 and 28. After comparison with A-DInSAR results and field inspection, we conclude that these are produced by very local movements related to monument instabilities [5].
The Advanced Differential Satellite Synthetic Aperture Radar (A-DInSAR) processing of the Lorca area Sentinel-1A images (see [5] for a detailed description of the advanced processing of the satellite radar images). Both orbits, ascending and descending (tracks 103 and 8 respectively) from Sentinel-1A, are used to decompose the measured LOS movement/mean velocities into horizontal (E-W) and vertical components [119][120][121] over the studied area. The A-DInSAR study covers the same time interval spanned by the GNSS campaigns (November 2015-February 2017). Radar data are processed using the Coherent Pixel Technique (CPT) [5,122]. The total area covered by the GNSS network is approximately 70 km 2 . A-DInSAR was processed for an extended region with a total area of 170 km 2 . In both geometries, ascending and descending, the study area is covered by three bursts of the same swath. A total of 42 Interferometric Wide Swath (IW) SLC images from the Sentinel-1A satellite resulted in 185 interferograms (22 images and 137 interferograms for ascending data, 20 and 48 respectively for descending [5]). The results are shown in Figure 23. Some selected time series are shown in the Supplementary Material by [5] for descending LOS. This is the first DInSAR study of the Alto Guadalentin Basin using two different geometries for the same time period.
Fernández et al. [5] obtain the vertical and E-W components of the displacement from the ascending and descending LOS motions. E-W and vertical components of the displacement fields in the area obtained using ascending and descending LOS results are shown in Figure 24. For this case [5], the decomposition above does not produce particularly good results due to methodological aspects and to the fact that the magnitude of the E-W motion is at sub-centimeter levels in many of the coherent pixels, i.e., the same order of the A-DInSAR uncertainty. The comparison of GNSS and A-DInSAR results allows to estimate the error of the A-DInSAR processing, relative to GNSS, as~0.7-1 cm [5].

Discussion and Inversion of the 3D Displacement Field
As previously noted, the establishment and observation (spanning November 2015-February 2017) of a local GNSS network allowed, for the first time, for measurement of the 3D displacement field in the Alto Guadaletin area, associated with exploitation of the local aquifer [5]. Also, for the first time, A-DInSAR results have been obtained using both ascending and descending radar images from the Sentinel-1A Copernicus radar satellite, allowing estimation of both vertical and horizontal (E-W) displacement components, a 2D displacement field, at higher spatial resolution than GNSS [5]. See Figures 23 and 24.
These results [5] highlight how the ad hoc establishment of survey mode GNSS networks improves the spatio-temporal monitoring of the 3D displacement field of areas subjected to extensive groundwater extraction, therefore representing a valuable monitoring technique. Moreover, GNSS observation provides complementary information to A-DInSAR results, allowing for their validation and scaling. In addition, at a local level, it is observed that the GNSS network does not completely cover the current displacement area, in particular along the SW region (Figures 23 and 24), because the network was defined based on displacements obtained prior to 2012 ( Figure 22). However, their Sentinel-1 A-DInSAR results [5] show that the deformation has extended in the SW direction, which today is the region of the most significant water extraction [123]. Therefore, the GNSS network needs to be extended over that region with additional GNSS stations.
The results obtained for the studied area highlight that [5]: (i) Simultaneous GNSS and A-DInSAR results are consistent with each other (Figure 24 and Supplementary Information from [5]; (ii) the results obtained for rates and pattern of the displacement are consistent with previous DInSAR results [21,22]; however, (iii), the horizontal displacement rate has a maximum amplitude of 2-3 cm/year ( Figure 23) and it is a significant component of the observed deformation field. Therefore, the horizontal displacement cannot be neglected, as it was done in previous studies [21,22].
Since the results by Fernández et al. [5] demonstrate that the horizontal displacements represent a significant component of the deformation field of the studied area, they carry out sensitivity tests by neglecting/including this horizontal motion in order to assess the variability (or bias percentage) on the determination of the aquifer characteristics and their temporal evolution using deformation modeling. To do this, they employ both the GNSS and A-DInSAR results. In addition, taking into account the linear time behavior of the displacement field [4], they also consider displacement rates in their study. Fernández et al. [5] use ten different data sets of surface displacement (denoted as Cases) covering the period November 2015 to February 2017 and carry out the inversion using the previously described forward model and inversion methodology. The cases are the following: (A) LOS A-DInSAR results obtained for descending orbit images, assuming 100% as vertical displacement; (B) LOS A-DInSAR results obtained for ascending orbit images, assuming 100% as vertical displacement; (C) Up-Down component obtained from the A-DInSAR, combining ascending and descending orbit images results; (D) Purely LOS A-DInSAR results obtained for descending orbit images; (E) Purely LOS A-DInSAR results obtained for ascending orbit images; (F) Purely LOS A-DInSAR results obtained for ascending and descending orbit images; (G) Up-Down and E-W A-DInSAR results obtained by combining ascending and descending orbit images results; (H) 3D displacements determined using the GNSS surveys results; (I) Purely LOS A-DInSAR results obtained for ascending and descending orbit images together with the 3D displacements determined using the GNSS surveys results; (J) Up-Down and E-W A-DInSAR results obtained by combining ascending and descending orbit images, together with 3D displacements determined using the GNSS surveys results.
Cases A-C are one-dimensional (1D), D-G are 2D (either indirectly by combining Up-Down and E-W in any of the two measured LOS, or directly supplying the values for both displacement components separately), H is a purely 3D dataset, and I and J are a combination of 2D and 3D data (2D+3D data).
A very important difference between the purely GNSS 3D data and the rest of these cases is the number of displacement rate data points [5]. For the 3D Case (H) data set, there are just 108 displacement rates, while for the remainder we have thousands of measurements (Table 3). Table 3. Numerical summary of the inversion results for the ten cases considered 1 [5]. Fernández et al. [5] inverted each case and estimate the volume changes of the water table (volume and geometry) assuming a given pressure change value. Moreover, based on hydrogeological observations, they impose the criteria that sources were shallower than one kilometer. A summary of the results obtained is provided in Table 3 and Figure 25. In these figures, the blue colors indicate negative pressure values, while white colors indicate positive pressure changes cells.

CASE
Inversion results include the volume and geometry of the active part of the aquifer which produced the measured displacements. Here this is quantified by the intensity, equal to the product of volume by pressure change, as it is impossible to determine both quantities separately. If we increase pressure, we decrease volume and vice versa. In order to determine a general geometry, they constrained the value of pressure change [1,17]. After a trial analysis, they considered a pressure value of -3 MPa [5], selecting the value that gives us a source geometry most consistent with the characteristics of the Lorca aquifer.
The inversion results obtained from these data sets can be organized into two subsets: (i) Cases A-C (1D, vertical displacement) and (ii) Cases D-J (2D and 2D+3D). Both sets of results are internally consistent, with respective scattering on the order of 9 and 6% respectively ( Table 3).
The results of (i) are, on average,~26% greater in intensity/volume than those of (ii), indicating that using only one component of the displacement field and assuming that the displacements are only vertical significantly overestimates the volume of water extracted during the study period (on the order of tens of hm 3 ). This can have an important effect in predictions of future volume variations and surface displacements.
The results for group (i) show the largest discrepancies between Case C and the other two cases (A and B). For Case C we obtain an intensity value~20% greater than the other ones. This data set has been obtained from combining LOS results from ascending and descending orbits and they suffer from the small systematic errors in the methodology [5]. This problem is reflected in two aspects, the poorer misfit values and the appearance of additional sources at the edges of the study area that try to adjust these errors ( Figure 25). As a result, it was concluded [5] that it is clearly more appropriate, in agreement with previous works [42], to mitigate this problem by using the ascending and/or descending LOS directly in the inversion procedure. As a result, Case C is not considered further. Cases A and B obtain very consistent results (Table 3), with a dispersion of~4% for group (i) and a mean intensity value of 40 MPa·km 3 [5].
For the group (ii) results, the data set for the vertical and E-W components determined using the combination of ascending and descending LOS (Cases G and J) ( Table 3) again provides worst results (higher intensities, more additional sources and poorer misfit values, Figure 25 and Table 3). Case H (only GNSS determined displacements) also provides higher misfits (Table 3) and poorer adjusted geometries (Figure 25h), probably a result of using data with several important limitations [5]. The first one is the reduced number of data points available in relation to the number of unknown parameters (approximately one hundred versus several thousands ( Table 3). As previously mentioned, in cases like this, a large number of measurement points are needed to reliably estimate the distribution of reservoir volume changes [5,48]. Also, their limited number produces a poorer spatial distribution than A-DInSAR data and does not cover the entire deformed area ( Figure 23). Finally, as described previously, in this reduced data set we have some anomalous results associated to very local effects (Figure 23b and Supplementary Information from [5]). To establish a denser and more extended GNSS network requires additional cost and time for field observation. A very important difference between the purely GNSS 3D data and the rest of these cases is the number of displacement rate data points [5]. For the 3D Case (H) data set, there are just 108 displacement rates, while for the remainder we have thousands of measurements (Table 3).
Fernández et al. [5] inverted each case and estimate the volume changes of the water table (volume and geometry) assuming a given pressure change value. Moreover, based on hydrogeological observations, they impose the criteria that sources were shallower than one kilometer. A summary of the results obtained is provided in Table 3     These positive pressure sources adjust the errors and the effects of other deformation sources, different from water extraction (e.g., of tectonic origin) [5].
Another important result [5] is that significant differences (3-4%, at the level of error or lower, see Table 3) are not observed between using just one LOS (ascending or descending), both LOS (ascending and descending) displacement data sets together, or both combined with the GNSS results. The estimation of source characteristics is very similar for all cases, just slightly changing the misfit of the minimum values for Cases D-F, I.
In summary, contrary to previous studies in the Lorca area, Fernández et al. [5] conclude that the measurement and use of the horizontal and vertical displacements at the surface is important for the prediction of future volume variations and surface displacements. These differences can have important effects on the design of monitoring systems, help in the decision-making process related to the sustainable management of the aquifer resources, and improve the assessment of potential hazards related to the aquifer exploitation. These main conclusions did not include any local assumptions which could condition our interpretation methodology [5], but they have general applicability worldwide. The Lorca case can be considered an extreme case, taking into account that it has significant E-W horizontal deformation, but only in geographically limited areas of maximum deformation. In other regions, where significant horizontal E-W deformation may be more scattered and cover more extended areas [5] the effect of considering vertical deformation alone could be even more dramatic.
In summary, Fernández et al [5] show, using inversion results from different remote sensing data sets and the described methodology for inversion, that the operational monitoring of the aquifer can be done using A-DInSAR with ascending and/or descending satellite radar images. GNSS, using continuous or survey observation mode, can be used for validation and scaling purposes. GNSS stations should be installed in those locations that will best constrain the A-DInSAR results in areas of zero deformation for reference, maximum deformation areas for scaling and study of the time variation, or in low coherence areas so that both techniques complement each other. Continuous GNSS observations are preferable, if possible. Their proposed methodology would potentially reduce the cost of the geodetic monitoring system in a very important way. This effective A-DInSAR monitoring can be accomplished with the freely-available Copernicus Sentinel-1A and -1B satellite data, considering their global coverage and repeatability, ensuring their effective use for monitoring on a global scale [5].

General Summary
We have described a new inversion methodology, first presented by Camacho et al. [1], for modeling space and terrestrial geodetic displacement and gravity changes data in active areas with application tests. We assume that surface deformation and gravity changes are due to density and/or pressure changes allocated within extended source structures below the surface.
The characteristics of this approach can be summarized as follows: (a) It permits simultaneous inversion of the 1D to 3D displacement data coming from different techniques, including terrestrial and space data (GNSS, DInSAR, leveling, etc.); (b) Non-planar, non-gridded, inexact data can be used; (c) It allows for objective modelling of two causative structures: Pressure and mass anomalies bodies; (d) Simple analytical and well-known expressions are used for direct calculation for mass and pressure cells and their aggregation; (e) This methodology works in a fully 3D context; (f) There is not additional a priori requirements on the geometry and types of the causative sources; (g) The method automatically determines the type and number, modeling the geometry, of the causative source structures; (h) A free geometry of the causative structures is described by aggregation of small elemental cells.
Surface deformation can be computed by adding the influence of the small considered prisms. Given that the used direct models are linear and the entire subsurface is assumed to be isotropic, superposition is allowable. Using this assumption, these linear equations permit the computation of surface deformation based on the superposition of many prismatic blocks within a compacting reservoir of any geometric shape [5]. The effect of each cell is computed using point sources [1,5].
When considering the validity of the results depending on the combination of the detected sources (if we have more than one pressure source), we must take into account the pressure change values employed and the adjusted distances between the different 3D pressure sources, which should be greater than 8 time their radii [124]. Solutions diverging from this criteria must be considering carefully and as a first, rough approach to the reality. Figure 26 show a schematic flow diagram of proposed inversion methodology.

3D grid of (empty) cells (point sources)
Gravity data dg + Covariance matrix Q g Growth process (scale factor) Extended 3D models for structures of anomalous density and/or pressure, defined as aggregation of (filled) cell/point sources Nearly automatic inversion process METHODOLOGY 3-D Geometry Figure 26. Schematic flow diagram of the described inversion methodology [1,5,17,53,54] where we start from the data set, the medium characteristics and its 3D gridding, to get (using the direct model equations and complementary conditions) a 3D source model of the anomalous sources via a growth process.

Conclusions
In active areas, the source modelling for geodetic (terrestrial and space) data fitting is usually based on a priori fixed geometries (spherical, ellipsoidal prolate or oblate, etc.) for point or extended sources/causative bodies characterized by anomalous density or pressure. This new approach avoids this limitation to model the volume of these anomalous extended bodies in a 3D context with a free geometry.
In a first version of the methodology, the Mogi model was used as analytical solution for the nucleus of strain model. The Mogi model is a linear, first order solution. In a second version of the methodology, Geertsma's nucleus of strain was implemented as the direct model solution, considering a poroelastic medium in place of a purely elastic one. This new version of the methodology allows us to consider the case of deformation associated with groundwater aquifer exploitation, increasing the applicability of the model/inversion methodology.
The described methodology is being implemented, PAF software [54], in the framework of the ESFRI (European Strategy Forum on Research Infrastructures) infrastructure EPOS [125] in the 3D-def service and will be available for a general use shortly.
Most of the active regions are characterized by complicated pattern of ground deformation resulting from multiple natural (e.g., volcanic loading producing inflation, deflation, dike intrusion, active faulting, flank instability and landslides) and anthropogenic (e.g., local uplift or subsidence) sources. We are developing an improvement of the described methodology which simultaneously estimates the typical ground deformation components in addition to the pressure source effects, allowing for a more general inverse methodology that includes dislocation sources, such as those given by the Okada's expressions [126].
Author Contributions: J.F. and A.G.C. have contributed equally in designing and writing this manuscript.