For example, in References [10, 11] the Ensemble Smoother (ES), that is, a Monte Carlo (MC)-based method, is used to estimate the parameters of a transversely isotropic and isotropic constitutive law by the joint assimilation of horizontal and vertical displacements. Later developments of the ES algorithm in the form of Multiple Data Assimilation (ES-MDA) are used in Reference [12] to constrain the reservoir compaction coefficient and the subsurface basement elastic modulus by assimilation of ascending and descending line-of-sight displacements revealed by Permanent Scatterer Interferometry (PSI). Recently, Reference [13] compares the results obtained from applying the ES and the ES-MDA to estimate the parameters of the modified Cam-Clay and the Vermeer-Neher (VN) visco-elasto-plastic models [3]. The growing interest on the ES as a data assimilation technique is mainly due to its straightforward implementation and the capability of using simultaneously the measurements collected at different times and locations in a 4-D assimilation. However, its convergence is guaranteed only with a large number of MC simulations that may be computationally unfeasible in demanding large-scale problems. Moreover, its effectiveness in constraining the model for a strongly non-linear dependence of the outcome on the uncertain parameters is theoretically questionable. This limitation can be overcome by employing the ES in a multiple assimilation fashion, with the introduction of the ES-MDA technique [14], but this implies a further significant increase of the computational cost for the overall procedure, with the ES algorithm repeatedly applied to different ensembles. It follows the need for fast model predictions that can be achieved by approximating the forward model operators with surrogate solutions, obtained at a highly reduced computational cost. In geomechanical applications, surrogate models have been recently employed for uncertainty quantification (UQ) and sensitivity analysis (SA). A functional data analysis technique is used to perform a vertical displacements UQ on a synthetic field-scale test by setting the oedometric compressibility as uncertain parameter [15]. Polynomial chaos proxies of the linear poroelastic problem with hydromechanical coupling allows a variance-based SA on Lame's constants, the Biot-Willis coefficient and the hydraulic mobility [16]. Another non-intrusive approach is used in sequentially coupled reservoir-geomechanical simulations to estimate the caprock safety factor in a probabilistic setting [17]. Global SA and parameter identification are carried out by means of Gaussian processes in a carbon capture and storage test case to model fault poromechanics and induced seismicity in a stochastic framework [18]. A similar approach is used in Reference [19] to investigate the onset of fault activation in deep hydrocarbon reservoirs with surrogates based on polynomial chaos expansion. Global SA are also employed to investigate geomechanical fractured reservoirs [20] and hydraulically fractured wells with reduced order models [21]. In this work, we focus our attention on the generalized polynomial chaos expansion (gPCE) [22, 23, 24, 25, 26]. Using a spectral representation in the random space, gPCE-stochastic solutions are expressed by orthogonal polynomials of the input parameters [27]. This approach has already proved successful in several different engineering applications [28, 29] with a growing interest in the Earth science community for UQ and SA [30, 31, 32, 33, 34, 35]. Moreover, gPCE is suitable for identification problems in the Bayesian inference framework as documented in previous studies [36, 37, 38, 39, 40], to cite just a few. Here we present an approach where we further extend the use of the gPCE as surrogate for regional geomechanical models with the aim of land subsidence prediction due to the exploitation of subsurface resources.