Giuseppe Gambolati

Mathematical Simulation of the Subsidence of Ravenna

Land subsidence at Ravenna is the result of aquitard and reservoir compaction caused, respectively, by extensive groundwater withdrawals from the unconsolidated Quaternary basin and gas production from a number of pre-Quaternary pools scattered over the area. Water pumpage paralleled the postwar industrial development of Ravenna until the middle seventies when consumption was drastically curtailed owing to the economic crisis and the activation of a new aqueduct. Gas production started in 1952. The exploitation of several reservoirs is currently under way and the search for new fields is still in progress. Geodetic records indicate that the maximum cumulative subsidence over the period 1950–1986, including a natural geologic settlement of perhaps 2 mm/yr, has been 1.30 m in the industrial zone of Ravenna. In 1980 the municipality promoted a reconnaissance study with the primary aim of providing the information base needed to reconstruct the actual occurrence, understand correctly the physical behavior and produce the essential input data to a mathematical model which realistically relates the subsidence of the city to groundwater withdrawal and gas removal with an emphasis on their respective influences. The results from the three-dimensional numerical simulations, performed with the aid of mixed finite element, finite difference and integral models, show that the primary responsibility for the regional land sinking should be placed on the subsurface water overdraft which occurred until the middle 1970s. Gas withdrawal plays a role restricted to the area overlying each reservoir with a magnitude depending on the depth of burial, thickness of mineralized rocks and overall volumetric production. A major environmental impact may be expected where the gas subsidence bowl is intersected by the Adriatic coastline.

Coupled and partially coupled Eulerian-Lagrangian Model of freshwater-seawater mixing

The problem of density-dependent transport of salt in unconfined coastal aquifers is solved numerically by means of an implicit Eulerian-Lagrangian finite element formulation. Such a formulation leads to symmetric positive definite finite element matrices which are ideally suited for efficient solution by preconditioned conjugate gradient methods. Additional known advantages of the formulation are unconditional stability, reduced numerical dispersion and suitability for parallel computation. The method has been used to study the effect of dewatering on seawater intrusion within a vertical cross section through an aquifer in southern Italy, related to the construction of a thermoelectric power plant. To investigate the extent to which the dependence of fluid density on salt concentration affects the numerical solution, the flow and advection-dispersion equations were solved in coupled (iterative), partially coupled (noniterative) and completely decoupled modes. Partial coupling was found to yield results very close to those obtained by full coupling but at great savings in computer time; the less rigorous decoupled approach led to results substantially different from those obtained through coupling and partial coupling. Effects of aquifer heterogeneity and the construction of a cutoff wall on seawater intrusion are discussed.

Importance of poroelastic coupling in dynamically active aquifers of the Po river basin, Italy.

Uncoupling between the flow field and the stress field in pumped aquifers is the basis of the classical groundwater hydrology. Recently, some authors have disputed the assumption of uncoupling with regard to both fluid dynamics and porous medium deformation. The issue is very important as it could undermine the traditional approach to simulate subsurface flow, analyze pumping tests, and predict land subsidence caused by fluid withdrawal. The present paper addresses the problem of coupling versus uncoupling in the Po river plain, a normally consolidated and normally pressurized basin which has experienced in the last 50 years a pronounced pore pressure drawdown because of water and gas removal and where a large hydromechanical database is available from the ground surface down to 4000 m depth. A numerical study is performed which shows that the matrix which relates flow to stress is very similar to the capacity matrix of the uncoupled flow equation. A comparison of results obtained with the finite element integration of the coupled and uncoupled models indicates that pore pressure is rather insensitive to coupling anywhere within the pumped formation while in the adjacent aquitard-aquifer units, coupling induces a slight overpressure which quickly dissipates in time with a small initial influence on medium deformation, and specifically on land subsidence. As a major consequence the uncoupled solutions to the fluid dynamic and the structural problems appear to be fully warranted on any timescale of practical interest in a typical normally consolidated and pressurized basin.

A fully coupled 3-D mixed finite element model of Biot consolidation

The numerical solution to the Biot equations of 3-D consolidation is still a challenging task because of the ill-conditioning of the resulting algebraic system and the instabilities that may affect the pore pressure solution. Recently new approaches have been advanced based on mixed formulations. In the present paper a fully coupled 3-D mixed finite element model is developed with the aim at alleviating the pore pressure numerical oscillations at the interface between materials with different permeabilities. A solution algorithm is implemented that takes advantage of the block structure of the discretized problem. The proposed model is verified against well-known analytical solutions and successfully experimented with in realistic applications of soil consolidation.

Ill-conditioning of finite element poroelasticity equations

Abstract The solution to Biots coupled consolidation theory is usually addressed by the finite element (FE) method thus obtaining a system of first-order differential equations which is integrated by the use of an appropriate time marching scheme. For small values of the time step the resulting linear system may be severely ill-conditioned and hence the solution can prove quite difficult to achieve. Under such conditions efficient and robust projection solvers based on Krylovs subspaces which are usually recommended for non-symmetric large size problems can exhibit a very slow convergence rate or even fail. The present paper investigates the correlation between the ill-conditioning of FE poroelasticity equations and the time integration step Δt. An empirical relation is provided for a lower bound Δtcrit of Δt below which ill-conditioning may suddenly occur. The critical time step is larger for soft and low permeable porous media discretized on coarser grids. A limiting value for the rock stiffness is found such that for stiffer systems there is no ill-conditioning irrespective of Δt however small, as is also shown by several numerical examples. Finally, the definition of a different Δtcrit as suggested by other authors is reviewed and discussed.

Is a simple diagonal scaling the best preconditioner for conjugate gradients on supercomputers

Abstract The implementation of accelerated conjugated gradients for the solution of large sparse systems of linear equations on vector/parallel processors requires programming features which are significantly different from the one needed on a scalar machine. Furthermore, a numerical algorithm which works well on the latter may be largely inefficient on the former. In the present analysis the numerical performances on a CRAY X-MP/48 of some widely known preconditioning techniques are compared, including the diagonal scaling, the incomplete Cholesky decomposition and the least squares polynomial preconditioners. The last ones, which are not suited to scalar machines, appear to be particularly attractive from a conceptual view point on vector/parallel architectures. The results obtained with 12 arbitrarily sparse finite element problems of growing size up to almost 5000, shows surprisingly that simple diagonal scaling is the most efficient preconditioning scheme in the majority of applications. In the few cases where it is not so, its performance is nevertheless comparable with that of the incomplete Cholesky factorization. Contrary to the general expectation, the polynomials exhibit a poor performance which does not increase with the degree and appear never to be competitive with the other two more traditional preconditioners.

A Block FSAI-ILU Parallel Preconditioner for Symmetric Positive Definite Linear Systems

A novel parallel preconditioner for symmetric positive definite matrices is developed coupling a generalized factored sparse approximate inverse (FSAI) with an incomplete LU (ILU) factorization. The generalized FSAI, called block FSAI, is derived by requiring the preconditioned matrix to resemble a block-diagonal matrix in the sense of the minimal Frobenius norm. An incomplete block Jacobi algorithm is then effectively used to accelerate the convergence of a Krylov subspace method. The block FSAI-ILU preconditioner proves superior to both FSAI and the incomplete block Jacobi by themselves in a number of realistic finite element test cases and is fully scalable for a given number of blocks.

Coastline regression of the Romagna Region, Italy, due to natural and anthropogenic land subsidence and sea level rise

The Romagna coastal area in the Northern Adriatic Sea has experienced in recent times continuous changes because of its precarious environment and low ground elevation above mean sea level (msl). Major processes that may influence the stability of the coast profile include land subsidence of both natural and anthropogenic origin and the msl rise caused by global climate change. According to the most accredited modeling predictions msl is expected to rise by almost 0.5 m over the next century because of the greenhouse effect. Natural land subsidence is the result of deep downward tectonic movement and consolidation of geologically recent deposits. It may be estimated in the range of 2–2.5 mm/yr in the Ravenna area and twice as much in the Po River delta. Anthropogenic land subsidence is primarily related to groundwater pumping from the upper fresh water aquifer system and gas production from Plio-Pleistocene reservoirs. Geodetic surveys from 1953 to 1990 provide documentary evidence of cumulative land settlement exceeding 0.8 m and 1.2 m at Marina di Ravenna and Cesenatico, respectively. In this study we estimate both natural and anthropogenic land subsidence for the years 2015, 2050, and 2100 with the aid of ad hoc finite element simulation models. The use of these predictions together with the expected msl rise shows that many present lowlands may be permanently submerged at the end of the next century. The extent of the flooded area of the Romagna coastal region can be as much as 690 and 910 km2, using optimistic and pessimistic land subsidence scenarios, respectively. A local detailed analysis indicates that the areas around the cities of Ravenna and Cesenatico may be seriously affected by sea water ingression while the city of Rimini is well protected because of its relatively high elevation above msl.

On the use of a main trend for the kriging technique in hydrology

Abstract A technique of interpolation based on a stochastic approach and referred to as ‘kriging’ technique has recently been contributed by the French School. A primary feature of the algorithm is its ability to provide an assessment of the predictive reliability. The accuracy of estimate depends on the evaluation of two stochastic quantities: the variogram γ and the main trend m of the hydrologic event z to be reconstructed. For an effective use of the method a correct understanding of the actual role played by m is required. With some ad hoc examples it is shown that using a polynomial trend with unspecified coefficients as suggested by the general theory may lead to paradoxical results whose behaviour is hard to predict a priori . It turns out that increasing the degree of m may yield an increase of the estimation error where one would expect to obtain the opposite. An alternative formulation is suggested which assumes m to be fully known in advance. Its expression is supposed to be derived from both the general behaviour of z as is recognizable from the available records and some extra-amount of information related to the general physical knowledge of the hydrological context. If this extra-amount of information is missing, the use of a constant trend should be recommended.

GIS simulations of the inundation risk in the coastal lowlands of the Northern Adriatic Sea

The Northern Adriatic Coastland, between the cities of Monfalcone and Cattolica, is characterized by locations of great tourist interest, such as the Venice Lagoon and the Romagna Riviera, and areas with a very precarious environmental setting, such as the Valli di Comacchio and the Po River Delta. Therefore, the coastal management and the design of new defence works of the littoral have to be made with the utmost care, possibly with the aid of numerical predictions of the coastal morphodynamics and the flood risk analysis of the lowland involved. In the study area, land may subside due to sediment natural compaction and subsurface fluid (water and gas) withdrawal. At the same time, littoral transport of solid material can contribute appreciably to change the shore morphology. Mean sea level may rise permanently due to global climate change (eustatism) and occasionally due to tides and intensive storm events. The predictions of each individual process is obtained using various ad hoc mathematical models and the outcome of the numerical simulations are managed with a GIS (geographical information system). Coastline evolution until the year 2100 is investigated and risk factor maps of the low-lying coastal areas are generated which account for the hazard of the expected event, and the land economic value and vulnerability.