On the Characterization and Forecasting of Ground Displacements of Ocean-Reclaimed Lands

: In this work, we study ground deformation of ocean-reclaimed platforms as retrieved from interferometric synthetic aperture radar (InSAR) analyses. We investigate, in particular, the suitability and accuracy of some time-dependent models used to characterize and foresee the present and future evolution of ground deformation of the coastal lands. Previous investigations, carried out by the authors of this paper and other scholars, related to the zone of the ocean-reclaimed lands of Shanghai, have already shown that ocean-reclaimed lands are subject to subside (i.e., the ground is subject to settling down due to soil consolidation and compression), and the temporal evolution of that deformation follows a certain predictable model. Speciﬁcally, two time-gapped SAR datasets composed of the images collected by the ENVISAT ASAR (ENV) from 2007 to 2010 and the COSMO-SkyMed (CSK) sensors, available from 2013 to 2016, were used to generate long-term ground displacement time-series using a proper time-dependent geotechnical model. In this work, we use a third SAR data set consisting of Radarsat-2 (RST-2) acquisitions collected from 2012 to 2016 to further corroborate the validity of that model. As a result, we veriﬁed with the new RST-2 data, partially covering the gap between the ENV and CSK acquisitions, that the adopted model ﬁts the data and that the model is suitable to perform future projections. Furthermore, we extended these analyses to the area of Pearl River Delta (PRD) and the city of Shenzhen, China. Our study aims to investigate the suitability of di ﬀ erent time-dependent ground deformation models relying on the di ﬀ erent geophysical conditions in the two areas of Shanghai and Shenzhen, China. To this aim, three sets of SAR data, collected by the ENV platform (from both ascending and descending orbits) and the Sentinel-1A (S1A) sensor (on ascending orbits), were used to obtain the ground displacement time-series of the Shenzhen city and its surrounding region. Multi-orbit InSAR data products were also combined to discriminate the up–down (subsidence) ground deformation time-series of the coherent points, which are then used to estimate the parameters ocean-reclaimed lands. In particular, our investigation enhances and extends the outcomes of other previous investigations to another area of interest, namely Shenzhen city, by studying the accuracy of the adopted ground deformation models in a more general context. A comparative analysis involving di ﬀ erent sets of independent SAR images has been performed, and the potential and limits of di ﬀ erent ground deformation models have been highlighted. We have shed light on the usefulness and applicability of the adopted models to forecast the future evolution of the terrain ground displacement. The results indicate that the geotechnical model used to characterize the deformation phenomena in the Shanghai area can also be successfully and in performed by other The results show that are more reasonable in prediction of deformation in Shenzhen. Taking the results of Shanghai into account, they allow us to say that the geotechnical model is applicable the description of the deformation in the reclamation area.


Introduction
Many coastal cities have implemented land reclamation projects, which can effectively provide feasible conditions for the development of coastal cities [1][2][3][4][5]. However, land reclamation projects artificially change the landform and geological environment, possibly resulting in geological disasters. Moreover, natural disasters such as flooding [6], coastal erosion [7], and ground subsidence in coastal areas are further aggravated by the current global sea-level rise [8]. Ground subsidence is an inherent problem in reclamation areas, due to the engineering and physical properties of the reclaimed materials. Severe ground subsidence could lead to the damage of buildings, roads, and public infrastructures. Therefore, the long-term investigation of the ground displacements of the reclamation areas is indispensable to maintain the public safety and sustainability of the coastal regions.
The Shanghai metropolitan area reclaimed a total of 580 km 2 from 1985 to 2010 [9]. Shenzhen, as an special economic zone, located on the south coast of China, has also reclaimed about 100 km 2 of lands from the ocean since 1979 [3]. Based on the engineering and geological properties of the materials used for the reclamation projects, the ground subsidence can be distinguished into three stages. Among them, the secondary compression stage of the alluvial deposits is the indispensable observation stage, which can last for more than ten years after the end of reclamation projects. Interestingly, earth observation (EO) systems have monitored the ongoing surface displacement phenomena for dozens of years [10]. Therefore, synthetic aperture radar (SAR) datasets acquired by different satellites provide us with an excellent opportunity to recover the spatial-temporal evolution of ground deformation of the reclaimed lands.
In previous work, Pepe et al. [23] applied a geotechnical model to link ENVISAT ASAR (ENV) (2007-2010) and COSMO-SkyMed (CSK) (2013-2016) ground deformation time-series of the Shanghai coastal area, thus generating long-term ENV+CSK surface displacement time-series [23]. However, there were no real InSAR-driven deformation measurements over the time gap between 2010 and 2013, but only information derived from the used geotechnical model to characterize the ground deformation signals from 2010 to 2013. The availability of additional SAR data covering the time gap is beneficial to shed light on the correctness of the ENV+CSK ground deformation time-series obtained, and the usefulness of the applied geotechnical model, to foresee the evolution of the ground displacement in the future, constrained by the available InSAR data. In particular, a new set of SAR data collected by the Radarsat-2 (RST-2) instrument, which is partially time-overlapped with the available CSK dataset (from January 2012 to October 2016), has been collected and processed. It is worth highlighting that the geotechnical model was built based on the specific dredger fill of the reclamation area of Shanghai [24]. Nevertheless, the applicability of such a geotechnical model to other reclamation areas and the rationality of the geotechnical model compared with other deformation models needs to be thoroughly investigated.
In this study, we performed an MT-InSAR analysis based on the application of the well-known small baseline subset (SBAS) technique [11] to the sets of SAR data available in both the Shanghai and Shenzhen areas. In particular, we first compared the RST2-SBAS and CSK-SBAS ground deformation time-series between December 2013 and March 2016 to check the agreement between the results of the two independent SAR datasets, assuming the ground deformation is mainly vertical (subsidence) Remote Sens. 2020, 12, 2971 3 of 24 in these areas, as demonstrated in the literature [32,33]. Then, we investigated the correctness of the geotechnical model used in Shanghai by comparing the combined ENV+CSK ground deformation time-series obtained using that model (see [23] for additional details) with those obtained by the independent set of RST-2 data. The engineering geology conditions of the Lingang New City, Shanghai were taken into account to identify the districts of the megacity that are more prone to subside, and be at risk of inundation in the future. In particular, we extend the work initially provided in [1,2,23] by testing the feasibility of the adopted model describing the ground subsidence against other potential time-dependent models.
Furthermore, we applied and extended these analyses to the area of the Pearl River Delta (PRD) and Shenzhen, China. First, by using two sets of ascending/descending ENV SAR data acquired between 2007 and 2010, we obtained the ground deformation time-series of the vertical (Up-Down) ground deformations by the use of the multi-temporal InSAR minimum acceleration combination (MinA) [34] technique. Then, we used the discriminated vertical ground deformations (i.e., subsidence) to foresee the temporal evolution of the deformation of the ocean-reclaimed platforms as dictated by three independent deformation models, including the geotechnical, the logarithmic, and the exponential decay models. Finally, we used a set of recently acquired Sentinel-1A (S1A) SAR data (ascending data track) to check the validity of the foreseen models. The results show that the independent SAR data and the forecasting models of the ground displacement are helpful for the continuous monitoring of coastal environments and the achievement of global sustainable development goals. There is evidence of a good agreement between the mean deformation velocities foreseen by the models constrained by the ENV data ten years ago and those obtained in the last year by the S1A data. The three models used are described in more detail in the paper.
The novel issues discussed in this paper concern: (i) the exploitation of InSAR data to constrain different ground deformation models, (ii) the cross-comparison and validation of the models used, (iii) the extension of the used models in another geographical area, and (iv) the characterization of the recent ground deformations of Shenzhen, China.

Lingang New City of Shanghai
Lingang New City is located in the southeast corner of Shanghai, bounded to the south by Hangzhou Bay and to the east by the East China Sea (approximately longitudes from 121.8 • E to 122 • E and latitudes from 30.83 • N to 31 • N), as shown in Figure 1a. An abundant supply of sediment from the Yangtze River made it possible for Shanghai to reclaim new lands, as shown in Figure 1a. In particular, the reclamation procedure of Lingang New City started in 1994 and was almost completed in 2006 [35]. The reclamation project has reclaimed more than 130 km 2 of new land for Lingang, which accounts for about 42% of the total area in Lingang New City. In addition to the previous area formed by sediment deposition before 1994, Lingang New City covers a total area of 315 km 2  Hangzhou Bay and to the east by the East China Sea (approximately longitudes from 121.8°E to 122°E and latitudes from 30.83°N to 31°N), as shown in Figure 1a. An abundant supply of sediment from the Yangtze River made it possible for Shanghai to reclaim new lands, as shown in Figure 1a. In particular, the reclamation procedure of Lingang New City started in 1994 and was almost completed in 2006 [35]. The reclamation project has reclaimed more than 130 km 2 of new land for Lingang, which accounts for about 42% of the total area in Lingang New City. In addition to the previous area formed by sediment deposition before 1994, Lingang New City covers a total area of 315 km 2 in the southeast corner of Pudong New District of Shanghai.

Shenzhen
Shenzhen is located in the south of Guangdong province and the eastern side of Pearl River Delta (PRD), bounded to the south by Shenzhen Bay and Hong Kong, to the east by Daya Bay, and to the west by the Pearl River estuary and Dapeng Bay, as shown in Figure 2a. For the sake of urban development, about 100 km 2 of lands have been reclaimed in Shenzhen since 1979 [3,36]. The consolidation settlement of the reclamation area also brings potential public safety risks to the Shenzhen coastal area.  (ascending passes, VV polarization, with an average revisit time of about 35 days, with a sensor sidelooking and a satellite heading angle of about 23° and -12°, respectively). The second dataset consists  of 44 images, collected by RST-2 sensor, from 29 June 2012, to 22 October 2016 (descending passes,  VV polarization, with an average revisit time of about 24 days, with a sensor side-looking and a  satellite heading angle of about 25° and 195°, respectively). The third dataset consists of 61 images, collected by CSK sensor, from 7 December 2013, to 18 March 2016 (descending passes, HH polarization, with an average revisit time of about 24 days, with a sensor side-looking and a satellite heading angle of about 29° and 8°, respectively). The spatiotemporal coverages of the three datasets are shown in Figure 1.

Shenzhen
Shenzhen is located in the south of Guangdong province and the eastern side of Pearl River Delta (PRD), bounded to the south by Shenzhen Bay and Hong Kong, to the east by Daya Bay, and to the west by the Pearl River estuary and Dapeng Bay, as shown in Figure 2a. For the sake of urban development, about 100 km 2 of lands have been reclaimed in Shenzhen since 1979 [3,36]. The consolidation settlement of the reclamation area also brings potential public safety risks to the Shenzhen coastal area.
Three different SAR datasets have been used to retrieve ground deformation of Shenzhen. The first dataset consists of 18 images, collected by ENV, from 13 June 2007, to 28 April 2010 (ascending passes, VV polarization). The second dataset consists of 20 images, also collected by ENV, from 24 June 2007, to 4 April 2010 (descending passes, VV polarization). The third dataset consists of 31 images, collected by S1A sensor, from 6 June 2019, to 12 June 2020 (ascending passes, VV polarization, with an average revisit time of about 12 days, with a sensor side-looking and a satellite heading angle of about 39° and 348°, respectively). The spatiotemporal coverages of those three datasets are shown in Figure 2.

InSAR Algorithms for the Estimation of the Ground Displacement over Ocean-Reclaimed Lands
For our experiments, the small baseline subset (SBAS) technique [11] was applied for the generation of the line-of-sight (LOS)-projected ground displacement time-series and mean deformation velocity maps [11]. Here, we briefly describe the rationale of the SBAS technique. The algorithm is a widely used multi-temporal InSAR (MT-InSAR) approach that is based on the process of long-term sequences of multi-look interferograms. In particular, if we suppose a set of Q co- Three different SAR datasets have been used to retrieve ground deformation of Shenzhen. The first dataset consists of 18 images, collected by ENV, from 13 June 2007, to 28 April 2010 (ascending passes, VV polarization). The second dataset consists of 20 images, also collected by ENV, from 24 June 2007, to 4 April 2010 (descending passes, VV polarization). Tables S4 and S5 in Supplementary Material list the ENV ascending and descending dataset respectively. The third dataset consists of 31 images, collected by S1A sensor, from 6 June 2019, to 12 June 2020 (ascending passes, VV polarization, with an average revisit time of about 12 days, with a sensor side-looking and a satellite heading angle of about 39 • and 348 • , respectively). Table S6 in Supplementary Material lists the S1A dataset. The spatiotemporal coverages of those three datasets are shown in Figure 2.

InSAR Algorithms for the Estimation of the Ground Displacement over Ocean-Reclaimed Lands
For our experiments, the small baseline subset (SBAS) technique [11] was applied for the generation of the line-of-sight (LOS)-projected ground displacement time-series and mean deformation velocity maps [11]. Here, we briefly describe the rationale of the SBAS technique. The algorithm is a widely used multi-temporal InSAR (MT-InSAR) approach that is based on the process of long-term sequences of multi-look interferograms. In particular, if we suppose a set of Q co-registered SAR images, acquired at ordered times t 1 , t 2 , . . . , t Q T , the SBAS algorithm requires a proper selection of SAR acquisitions of data pairs characterized by short temporal (i.e., the time interval between two acquisitions) and spatial (i.e., the distance between two satellite orbits) baselines. The implementation of such constraints on the baselines allows one to reduce the noise related to decorrelation effects in the interferograms and, accordingly, to maximize the number of reliable measurement points.
On the other hand, such a preliminarily pre-selection step of the small baseline (SB) interferograms may lead to the SAR images being distributed into subsets, which are separated by large baselines and are independent each other, i.e., with no eligible interferograms connecting data belonging to different subsets. This circumstance gives rise to an undetermined problem for the ground deformation time-series inversion, which is addressed by searching for a minimum norm least-squares (LS) solution through the singular value decomposition (SVD) method [37]. In our study, we use an improved SBAS processing chain, which also includes the detection and correction of residual topographic errors related to the inaccurate digital elevation models (DEMs) used for computing the differential interferograms and the effects of orbital inaccuracies. The phase unwrapping (PhU) operations are performed with the well-known space-time extended minimum cost flow (EMCF) PhU technique initially presented in [38]. Finally, the obtained LOS-projected ground deformation time-series are cleaned from the atmospheric phase screen (APS) signals. Additional details on the SBAS algorithm and the SBAS processing chain, can be found in [39]. After the SBAS processing, the LOS ground displacement time-series related to different multiple-orbit (ascending/descending) data tracks are geocoded to a common spatial grid of high-coherent points. For our subsequent analyses, let d(P) = d 1 (P), d 2 (P), . . . , d Q (P) T be the (geocoded) LOS deformation time-series, which is computed in correspondence with the generic pixel P. Therefore, a post-processing step, aimed at properly combining them to retrieve the East-West (E-W), and North-South (N-S), and up-down (U-D) components of the surface displacement [34], is applied. Indeed, LOS ground deformation values can be expressed as [34] d(P) = sin ϑ(P) cos ϕ(P) · e(P) + sin ϑ(P) sin ϕ(P) · n(P) + cos ϑ(P) · h(P) (1) where e(P) = e 1 (P), e 2 (P), . . . , e Q (P) T , n(P) = n 1 (P), n 2 (P), . . . , n Q (P) T and h(P) = h 1 (P), h 2 (P), . . . , h Q (P) T are the 3-D cumulative deformation components in the East-West, North-South, and vertical directions, respectively. Moreover, ϑ(P) is the radar-to-target incidence angle and ϕ(P) is the corresponding satellite-heading angle. The geometry of the multiple-angle SAR acquisitions (ascending/descending) is shown in Figure 3. The adopted combination technique is known as the minimum acceleration (MinA) technique since it relies on imposing the constraint that the (unknown) 3-D deformation time-series has minimum acceleration. To decompose the preliminarily computed, LOS-projected, grounddisplacement time-series into the 3-D ground-displacement time-series, we adopted the procedure fully detailed in [34]. It consists of connecting, pixel-by-pixel, the LOS-projected displacement timeseries to the set of (unknown) 3-D velocities among consecutive SAR acquisitions, and the subsequent time integration of the 3-D velocities [34]. Mathematically, this operation translates into writing an underdetermined system of independent linear equations, which can be solved using LS minimization based on the SVD method. The obtained system of linear equations is regularized by imposing the additional constraint that the 3-D displacement components have minimum acceleration, i.e., that the difference between consecutive velocities is minimal. Finally, once estimated, the unknown 3-D velocities are time-integrated to recover the E-W, N-S, and U-D deformation time series. The capability of the novel MinA technique to allow detailed Earth surface deformation analyses is fully exploited here by profitably combining independent datasets acquired by different satellite SAR sensors (ENVISAT, RADARSAT-2, Cosmo-SkyMed) operating at different wavelengths (i.e., C, X), polarization (i.e., HH, VV), and looking angles [40]. In particular, in our paper, we applied a simplified version of the MinA technique, characterized by assuming the North-South components of the ground deformation are negligible in the measured LOS deformation timeseries. Interested readers are invited to refer to [1] for a fully detailed characterization of the modified MinA approach, where an extensive statistical study of the error budget of the method is also presented.

Ground Deformation Results of Shanghai and Analyses
This section is devoted to proving that the geotechnical model used to fill the time gap between the ENV and the CSK deformation time-series in [23] is sound. To this aim, the independent RST-2 SAR dataset and the ground deformation time-series, obtained by processing this new SAR dataset with the SBAS algorithm [11], represent the key to support this proof. As shown in Figure 4, the RST-2 and CSK datasets have a long overlap period from December 2013 to March 2016, and the RST-2 The adopted combination technique is known as the minimum acceleration (MinA) technique since it relies on imposing the constraint that the (unknown) 3-D deformation time-series has minimum acceleration. To decompose the preliminarily computed, LOS-projected, ground-displacement time-series into the 3-D ground-displacement time-series, we adopted the procedure fully detailed in [34]. It consists of connecting, pixel-by-pixel, the LOS-projected displacement time-series to the set of (unknown) 3-D velocities among consecutive SAR acquisitions, and the subsequent time integration of the 3-D velocities [34]. Mathematically, this operation translates into writing an underdetermined system of independent linear equations, which can be solved using LS minimization based on the SVD method. The obtained system of linear equations is regularized by imposing the additional constraint that the 3-D displacement components have minimum acceleration, i.e., that the difference between consecutive velocities is minimal. Finally, once estimated, the unknown 3-D velocities are time-integrated to recover the E-W, N-S, and U-D deformation time series. The capability of the novel MinA technique to allow detailed Earth surface deformation analyses is fully exploited here by profitably combining independent datasets acquired by different satellite SAR sensors (ENVISAT, RADARSAT-2, Cosmo-SkyMed) operating at different wavelengths (i.e., C, X), polarization (i.e., HH, VV), and looking angles [40]. In particular, in our paper, we applied a simplified version of the MinA technique, characterized by assuming the North-South components of the ground deformation are negligible in the measured LOS deformation time-series. Interested readers are invited to refer to [1] for a fully detailed characterization of the modified MinA approach, where an extensive statistical study of the error budget of the method is also presented.

Ground Deformation Results of Shanghai and Analyses
This section is devoted to proving that the geotechnical model used to fill the time gap between the ENV and the CSK deformation time-series in [23] is sound. To this aim, the independent RST-2 SAR dataset and the ground deformation time-series, obtained by processing this new SAR dataset with the SBAS algorithm [11], represent the key to support this proof. As shown in Figure 4, the RST-2 and CSK datasets have a long overlap period from December 2013 to March 2016, and the RST-2 dataset partially fills the gap between the ENV and the CSK datasets, from January 2012 to December 2013. To compare the independent CSK and RST-2 InSAR ground deformation time-series, we assumed that the deformation of the terrain in the coastal regions of the city is mainly vertical (e.g., subsidence), as also assumed in previous investigations [19] and supported by evidence in the literature [32]. Under this hypothesis, the vertical ground deformation was obtained by merely re-projecting the line-of-sight (LOS) deformation into the vertical direction, as clarified in the following sub-section.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 25 Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing dataset partially fills the gap between the ENV and the CSK datasets, from January 2012 to December 2013. To compare the independent CSK and RST-2 InSAR ground deformation time-series, we assumed that the deformation of the terrain in the coastal regions of the city is mainly vertical (e.g., subsidence), as also assumed in previous investigations [19] and supported by evidence in the literature [32]. Under this hypothesis, the vertical ground deformation was obtained by merely reprojecting the line-of-sight (LOS) deformation into the vertical direction, as clarified in the following sub-section.

RST-2 Ground Deformation Time-Series: As the Third Party Inspection Data
The RST-2 dataset was processed using the SBAS algorithm [11] to obtain the line-of-sight (LOS)projected displacement time-series, which were preliminarily geocoded to have a common spatial reference grid over which to perform the subsequent analyses. First, we generated 145 RST-2 differential SAR interferograms, selected by considering maximum perpendicular and temporal baseline values of 400 m and 180 days, respectively. Moreover, a complex multi-look operation, with four looks in the azimuth and twenty looks in the range direction, was performed for every interferogram to reduce the influence of the decorrelation noise. Correspondingly, we obtained 145 multi-looked interferograms with a resolution of about 90 × 90 m. All the interferograms were noisefiltered [41] and then unwrapped [38,42]. After that, we obtained the RST-2 LOS-projected displacement time-series. As earlier anticipated, it is well known that deformation is mainly vertical in Shanghai [32]. Thus, the RST-2 LOS displacement time-series were converted, pixel-by-pixel, into vertical movements as follows: where is the vertical deformation, is the LOS deformation, and is the radar side-looking angle. The CSK and RST-2 mean deformation velocity maps of the Lingang New City are shown in Figure 5.

RST-2 Ground Deformation Time-Series: As the Third Party Inspection Data
The RST-2 dataset was processed using the SBAS algorithm [11] to obtain the line-of-sight (LOS)-projected displacement time-series, which were preliminarily geocoded to have a common spatial reference grid over which to perform the subsequent analyses. First, we generated 145 RST-2 differential SAR interferograms, selected by considering maximum perpendicular and temporal baseline values of 400 m and 180 days, respectively. Moreover, a complex multi-look operation, with four looks in the azimuth and twenty looks in the range direction, was performed for every interferogram to reduce the influence of the decorrelation noise. Correspondingly, we obtained 145 multi-looked interferograms with a resolution of about 90 × 90 m. All the interferograms were noise-filtered [41] and then unwrapped [38,42]. After that, we obtained the RST-2 LOS-projected displacement time-series. As earlier anticipated, it is well known that deformation is mainly vertical in Shanghai [32]. Thus, the RST-2 LOS displacement time-series were converted, pixel-by-pixel, into vertical movements as follows: where h is the vertical deformation, d is the LOS deformation, and ϑ is the radar side-looking angle. The CSK and RST-2 mean deformation velocity maps of the Lingang New City are shown in Figure 5. Initially, we compared the RST-2 (vertical) surface deformation time-series with that provided by the (single) CSK displacement time-series during the overlapping time, i.e., from December 2013 to March 2016, to evaluate the consistency of the InSAR measurements obtained by the independently processed CSK and RST-2 SAR data. However, of course, the time acquisitions T1 ≡ [t(1), t(2), · · · , t(i)] (i = 1, 2, · · · , M) of the CSK data are not the same as those of the RST-2 dataset, namely T2 ≡ [t(1), t(2), · · · , t(i)] (i = 1, 2, · · · , N), with t(i), with the time of the i-th SAR image for the two separate datasets being, from December 2013 to March 2016. Accordingly, we preliminarily (linearly) time-interpolated the CSK deformation time-series over the RST-2 acquisition times as follows: where d CSK (t(i)) is the CSK cumulative ground deformation at t(i), where t(i) belongs to T2. Moreover, d CSK (tr(i)) and d CSK (tl(i)) are the CSK cumulative deformation at tr(i) and tl(i), respectively, with tr(i) and tl(i) being the closest times to the right and left of t(i), respectively. As a result, we obtained the interpolated CSK ground deformation time-series at the time acquisition vectors T2. However, the ground deformation time-series for every dataset is calculated relative to the earliest SAR image of the set (the earliest RST-2 image is 29 January 2012, and the earliest CSK image is 7 December 2013). Accordingly, to compare the two CSK and RST-2 datasets, the systematic ground deformation bias between the two SAR datasets was estimated for every point that was common to both datasets. Expressed in a mathematical formula, the following equation holds (see also Figure 4): where d CSK (P,t) is the CSK deformation time-series at the given point P into the overlapped spatial region between the two processed CSK and RST-2 slices (see Figure 1a), d ∼ RST−2 (P,t) is the adjusted RST-2 ground-deformation time-series in the same point P, and d RST−2 (P,t) is the corresponding, original RST-2 ground deformation time-series. Moreover, we calculated, for every point P, the unique value of the ground-deformation bias ∆ that minimizes d CSK − d ∼ RST−2 2 . The mean value of root mean square error (RMSE) between the adjusted CSK and RST-2 ground deformation time-series, calculated at common coherent points between d CSK and d ∼ RST2 , was estimated. Its average value over the scene is about 3.5 mm, which is in agreement with the expected accuracy of SBAS measurements, which is in the order of 5 mm for the ground-deformation time-series [43]. We also cross-compared the (vertical) mean ground deformation velocity values over the common coherent points of the RST-2 and CSK datasets, considering the overlapped period from December 2013 to March 2016. A consistency analysis between the two velocities maps was performed by least-square estimation using the total shared high-coherent pixels. The scatter plot is shown in Figure 6, which also portrays the corresponding mathematical expressions of the regression line between the RST-2 and the CSK ground-deformation velocities. Considering also the values of the ground deformation mean velocity, the agreement of the adjusted independent CSK and RST-2 SAR datasets is consistent.   March 2016. A consistency analysis between the two velocities maps was performed by least-square estimation using the total shared high-coherent pixels. The scatter plot is shown in Figure 6, which also portrays the corresponding mathematical expressions of the regression line between the RST-2 and the CSK ground-deformation velocities. Considering also the values of the ground deformation mean velocity, the agreement of the adjusted independent CSK and RST-2 SAR datasets is consistent.

Analysis of Time-Gapped ENV+CSK and RST-2 Ground Displacement Time-Series
In the literature [23], a method to combine the time-gapped ENV and CSK displacement timeseries, benefiting from the knowledge of a proper time-dependent geotechnical model, was initially developed. The method relied on the estimation of the model parameters by using the Levenberg-Marquadt optimization [44,45] to constrain the model to the ENV and CSK data simultaneously. In particular, the adopted geotechnical model was derived from laboratory tests that simulated, at a reduced scale, real scenarios and the kinematic parameters that could be related to dredger-fill soil characteristics (i.e., thickness, water content, void ratio). The time-dependent model has hyperbolic behavior, and its mathematical expression is as follows: where t is the vector of the ordered times of the available SAR acquisitions, S0 is the asymptotic value of vertical deformation (theoretically assumed at infinite time), k and are two parameters that control the shape of a given curve among the family of all possible curves, and is a time-delay, which takes into account the uncertainty of the knowledge of the exact time when reclamation processes ended. Figure 7 reproduces the ground mean deformation velocity map of the Lingang New City from February 2007 to March 2016, as obtained by using the model (5) and the method to

Analysis of Time-Gapped ENV+CSK and RST-2 Ground Displacement Time-Series
In the literature [23], a method to combine the time-gapped ENV and CSK displacement time-series, benefiting from the knowledge of a proper time-dependent geotechnical model, was initially developed. The method relied on the estimation of the model parameters by using the Levenberg-Marquadt optimization [44,45] to constrain the model to the ENV and CSK data simultaneously. In particular, the adopted geotechnical model was derived from laboratory tests that simulated, at a reduced scale, real scenarios and the kinematic parameters that could be related to dredger-fill soil characteristics (i.e., thickness, water content, void ratio). The time-dependent model has hyperbolic behavior, and its mathematical expression is as follows: where t is the vector of the ordered times of the available SAR acquisitions, S 0 is the asymptotic value of vertical deformation (theoretically assumed at infinite time), k and λ are two parameters that control the shape of a given curve among the family of all possible curves, and δ is a time-delay, which takes into account the uncertainty of the knowledge of the exact time when reclamation processes ended. Figure 7 reproduces the ground mean deformation velocity map of the Lingang New City from February 2007 to March 2016, as obtained by using the model (5) and the method to link time-gapped data presented in [23]. Note also that the model (5) represents a generalization of the model adopted in [46]. The RST-2 SAR data, however, partially time-overlap the CSK data sets (see the red region highlighted in Figure 4). This circumstance has led us to the idea of testing the validity of the retrieved best-fit models using this new independent set of RST-2 data. To this aim, we estimated and compensated, for every coherent point common to the ENV+CSK and RST-2 datasets, the ground deformation bias, namely Ω, existing between the best-fit model d model [23] and the new RST-2 ground deformation time-series in the period from January 2012 and December 2013. It is evident (see also Figure 4) that: whered RST−2 is the corrected RST-2 deformation time-series adjusted to the model-driven ENV+CSK ground-deformation time-series already obtained and discussed in [23]. The bias Ω was estimated, for every point common to both SAR datasets, using a simple least-squares (LS) adjustment of the following problem: Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing ( , ) = ( , ) − Ω( ) where ^ is the corrected RST-2 deformation time-series adjusted to the model-driven ENV+CSK ground-deformation time-series already obtained and discussed in [23]. The bias was estimated, for every point common to both SAR datasets, using a simple least-squares (LS) adjustment of the following problem: Then, we calculated the RMSE between the adjusted RST-2 ground deformation time-series and the ENV+CSK displacement time-series at each common coherent point, which is shown in Figure 8. The results demonstrate, except for some noisy areas in the coastal regions afflicted by severe temporal decorrelation artifacts, that the ENV+CSK and RST-2 ground deformation products are in good agreement, thus making it evident that the model used in [23] was suitably correct.  Then, we calculated the RMSE between the adjusted RST-2 ground deformation time-series and the ENV+CSK displacement time-series at each common coherent point, which is shown in Figure 8. The results demonstrate, except for some noisy areas in the coastal regions afflicted by severe temporal decorrelation artifacts, that the ENV+CSK and RST-2 ground deformation products are in good agreement, thus making it evident that the model used in [23] was suitably correct.

Engineering Geology Analyses
In this section, we analyze the ground deformation results derived by the SBAS and the geotechnical model, taking into account the engineering geology condition of Lingang New City. From 1994 to 2006, a series of reclamation projects were carried out from the west side to the east side of Lingang New City. Abundant dredger fill, which is a sort of unconsolidated soil with high compressibility, large void ratio, and high-water content, has been primarily filled in the intertidal areas along the Lingang New City coast.
Based on the geological environment and the magnitude of deformation velocity, Lingang New City can be roughly divided into three zones (A, B, and C), as shown in Figure 9. Zone A is located in the western part of the city, and it was reclaimed before 1973. The ground settlement of Zone A has been unremarkable. Zone B is located in the central part of Lingang New City, and it was reclaimed between 1973 and 1994. This area has experienced a phase of rapid consolidation. Zone C, which is the eastern part of Lingang New City, has been newly reclaimed since 2002. Unconsolidated dredger fill is primarily distributed in zone C. Our results show that Zone A is stable after about a half century consolidation time, Zone B has small ground subsidence values after two decades of the consolidation phase, and the newly reclaimed Zone C is still experiencing severe ground subsidence.

Engineering Geology Analyses
In this section, we analyze the ground deformation results derived by the SBAS and the geotechnical model, taking into account the engineering geology condition of Lingang New City. Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing has been unremarkable. Zone B is located in the central part of Lingang New City, and it was reclaimed between 1973 and 1994. This area has experienced a phase of rapid consolidation. Zone C, which is the eastern part of Lingang New City, has been newly reclaimed since 2002. Unconsolidated dredger fill is primarily distributed in zone C. Our results show that Zone A is stable after about a half century consolidation time, Zone B has small ground subsidence values after two decades of the consolidation phase, and the newly reclaimed Zone C is still experiencing severe ground subsidence. We also calculated the mean RMSE values between the modeled ENV+CSK and the adjusted RST-2 deformation time-series over the three zones. The RMSE values are listed in Table 1.
We also compared the ground measurements obtained by the geotechnical model and RST-2 SBAS-InSAR and the engineering geology condition on the profile from I to I', which was derived from eight boreholes. The boreholes have a depth range from 30 to 50 m (as shown in Figure 9) and were collected to show the geological features of soil layers. The sub-soil in Lingang New City can be divided into six engineering geological layers, based on the geological ages and sub-soil properties. We now discuss the map shown in Figure 10c. The dredger fill, layer ① 3, is unevenly distributed in Lingang New City. It is mainly distributed in Zone B with a thickness not exceeding 3 We also calculated the mean RMSE values between the modeled ENV+CSK and the adjusted RST-2 deformation time-series over the three zones. The RMSE values are listed in Table 1. We also compared the ground measurements obtained by the geotechnical model and RST-2 SBAS-InSAR and the engineering geology condition on the profile from I to I', which was derived from eight boreholes. The boreholes have a depth range from 30 to 50 m (as shown in Figure 9) and were collected to show the geological features of soil layers. The sub-soil in Lingang New City can be divided into six engineering geological layers, based on the geological ages and sub-soil properties.
We now discuss the map shown in Figure 10c. The dredger fill, layer 1 3, is unevenly distributed in Lingang New City. It is mainly distributed in Zone B with a thickness not exceeding 3 m and evenly distributed in Zone C with a thickness ranging from 3 to 6 m. The silty clay, layer 2 1, is relatively stable compared with 1 3. It is mainly distributed in the west and north of Lingang New City. The sandy silt of layer 2 3 has similar physical and mechanical properties as the dredger fill layer 1 3. It is prone to compression, and is distributed in the whole area of Lingang New City. In addition, the sandy silt of layer 2 3 is apt to form quicksand. The muddy clay layer 4 and the clay 5 1 have high compressibility and are representative soft soil layers in Shanghai. Layers 4 and 5 1 are sensitive to large-scale construction. Irregular construction may lead to an overwhelming ground deformation. However, the hard geological layers 6 and 7 are the stable soil layers. The consolidation of layers 1 3, 2 , 4 , and 5 could mainly result in ground deformation [22]. Remote Sens. 2020, 12,

Ground Deformation Results of Shenzhen and Analyses
For the monitoring of reclaimed land regions, some time-dependent models have been developed to follow the ongoing evolution and foresee the future evolution of the ground deformations. In addition to the geotechnical model verified in this study, we also explored the suitability of other two time-dependent models with i) logarithmic and ii) exponential decay [26,48] to be used for ground deformations monitoring and forecasting. However, it is first worth exploring whether the geotechnical model used to characterize the area of Shanghai is still valid for other landreclamation areas. To achieve our goal, we took Shenzhen, located east of the Pearl River Delta (PRD) We plotted in Figure 10 the absolute values of the mean ground deformation velocity obtained by the geotechnical model constrained by the time-gapped ENV and CSK SBAS-InSAR deformation time-series and the RST-2 SBAS-InSAR deformation measurements, as well as engineering geology of the profile from I to I'. Three boreholes, labeled as LGL30, LGL31, and LGL32, are located in Zone B. The other five boreholes, labeled as LGL33, LGC74, LGG70, LGC75, and LGL34, are located in the newly reclaimed Zone C. As shown in Figure 10, the velocity differences of the portion from LGL30 to LGL33 do not exceed 1 mm/year. Correspondingly, the thickness of the dredger fill layer in the portion from LGL30 to LGL33 does not exceed 2 m. The mean ground deformation velocity difference of the portion from LGL33 to LGC74 is smaller than 2 mm/year. This portion is the central part of the newly reclaimed area Zone C. As shown in Figure 10, along the profile the deformation velocity derived by both the geotechnical model and RST-2 SBAS-InSAR show an increasing trend. The portion from LGC74 to LGG75 has the thickest dredger fill. Correspondingly, it has relatively high deformation velocities and the highest ground deformation velocity differences. The deformation difference is higher than 3 mm/year and smaller than 5 mm/year. Since the portion from LGC75 to LGL34 is still under reclamation, although this portion has relatively thin dredger fills, it has the most significant ground subsidence values of the whole profile. The ground deformation difference values of the portion from LGC75 to LGL34 are smaller than 2 mm/year.

Ground Deformation Results of Shenzhen and Analyses
For the monitoring of reclaimed land regions, some time-dependent models have been developed to follow the ongoing evolution and foresee the future evolution of the ground deformations. In addition to the geotechnical model verified in this study, we also explored the suitability of other two time-dependent models with i) logarithmic and ii) exponential decay [26,48] to be used for ground deformations monitoring and forecasting. However, it is first worth exploring whether the geotechnical model used to characterize the area of Shanghai is still valid for other land-reclamation areas. To achieve our goal, we took Shenzhen, located east of the Pearl River Delta (PRD) region in China, into consideration. Since 1979, the coastal area of Shenzhen, which extends more than 100 km 2 (as shown in Figure 2a), has been modified in its geomorphology [3,36] due to land reclamation activities. Thus, they are anthropogenic at most [49,50]. In this section, we test the validity of such different reclamation models in the Shenzhen area. To this aim, three SAR datasets (the ascending ENV, the descending ENV, and the S1A) were used.

The Ground Deformation Results of Shenzhen City
The ascending ENV (A-ENV), descending ENV (D-ENV), and S1-A datasets were processed separately using the SBAS algorithm [11] to obtain the related line-of-sight (LOS)-projected deformation time-series. Accordingly, we have computed the mean ground deformation velocity maps for the A-ENV and D-ENV SAR datasets, which are shown in Figure 11. Subsequently, we have applied the multi-satellite MinA method [34] to recover the East-West (E-W) and the up-down (U-D) ground-deformation time-series. This procedure is fundamental for the subsequent analyses because it allows us to discriminate LOS-projected ground deformation measurements from the ground-subsidence signals (i.e., the U-D components of the terrain displacement) that are strictly related to the process of consolidation of the reclaimed lands. Figure 12 portrays the maps of 2-dimensional ground deformation velocity components. The areas with sensible vertical deformation (subsidence) are almost located in the reclamation regions. To correlate the identified deformation signals to the process of consolidation of the reclaimed terrains, we isolated and focused on the areas characterized by significant up-down ground deformations, from June 2007 to April 2010, as shown in Figure 13a. We masked those pixels whose absolute value of the East-West ground deformation was higher than 10 mm/year. Unfortunately, no descending S1A data are available for the Shenzhen area from June 2019 to June 2020. Thus, we can not retrieve the two-dimensional deformation of the terrain from June 2019 to June 2020 by the MinA method. To circumvent this problem, we have assumed that areas with no significant lateral movements in 2009 preserve this characteristic over time, and we have retrieved the vertical deformation of the terrain by simply re-projecting the LOS ground deformation of S1A using Equation (2) and the relevant side-looking angle (ϑ = 39 • ) of the Sentinel platform S1A sensor. As a result, Figure 13b portrays the map of S1A up-down ground deformation velocity components from June 2019 to June 2020.
Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing Sentinel platform S1A sensor. As a result, Figure 13b portrays the map of S1A up-down ground deformation velocity components from June 2019 to June 2020. We selected four points labelled as P1, P2, P3, and P4 to describe the ENV and S1A vertical deformation time-series, as shown in Figure 14. We can see that the deformation is relatively broad from June 2007 to April 2010. In contrast, the ground deformation slows down from June 2019 to June 2020. Figure 11. Shenzhen area results. Maps of ascending (left) and descending (right) mean deformation velocities retrieved by the ascending and descending ENV datasets. Only well-processed coherent SAR pixels exhibiting a temporal coherence larger than 0.6 were considered. Figure 11. Shenzhen area results. Maps of ascending (left) and descending (right) mean deformation velocities retrieved by the ascending and descending ENV datasets. Only well-processed coherent SAR pixels exhibiting a temporal coherence larger than 0.6 were considered.  We selected four points labelled as P1, P2, P3, and P4 to describe the ENV and S1A vertical deformation time-series, as shown in Figure 14. We can see that the deformation is relatively broad from June 2007 to April 2010. In contrast, the ground deformation slows down from June 2019 to June 2020.

Comparison and Verification of the Geotechnical Model and Other Deformation Models in Shenzhen Reclamation Area
It is well known from the literature [26,27,46,48,[51][52][53] that the consolidation process of a saturated clay layer in reclaimed lands involves two main stages: the primary consolidation and the secondary compression [48]. The primary consolidation is responsible for the most significant proportion of the total ground displacement that evolves quickly just after the end of reclamation projects. Conversely, the secondary compression and fill creep can last for decades after the end of the reclamation procedures and evolves slowly over time. This phenomenon can be explained by the well-known Terzaghi theory of consolidation [54] in terms of dissipation of excess pore pressures (i.e., the increase of significant stress) induced by the external loading in undrained conditions. Among others, the geotechnical model of Equation (5) was obtained by Yang et al. [24] from laboratory centrifuge tests based on the dredger fill characteristics of Shanghai. Through previous studies [19,23] and the verification analysis in Section 3 of our study, the model can commendably describe the secondary compression process of the dredger fill in the reclamation area of Shanghai. Moreover, the logarithmic model [26,48] and the exponential decay model [26] were used in the investigation of the long-term reclamation settlement of Hong Kong Airport. In the following, these two additional models are presented and discussed. Remote Sens. 2020, 12,

Comparison and Verification of the Geotechnical Model and Other Deformation Models in Shenzhen Reclamation Area
It is well known from the literature [26,27,46,48,[51][52][53] that the consolidation process of a saturated clay layer in reclaimed lands involves two main stages: the primary consolidation and the secondary compression [48]. The primary consolidation is responsible for the most significant proportion of the total ground displacement that evolves quickly just after the end of reclamation projects. Conversely, the secondary compression and fill creep can last for decades after the end of the reclamation procedures and evolves slowly over time. This phenomenon can be explained by the well-known Terzaghi theory of consolidation [54] in terms of dissipation of excess pore pressures (i.e., the increase of significant stress) induced by the external loading in undrained conditions. Among others, the geotechnical model of Equation (5) was obtained by Yang et al. [24] from laboratory centrifuge tests based on the dredger fill characteristics of Shanghai. Through previous studies [19,23] and the verification analysis in Section 3 of our study, the model can commendably describe the secondary compression process of the dredger fill in the reclamation area of Shanghai. For a constant secondary compression index, the secondary compression displacement [48] can be obtained as: where H 0 is the initial thickness of the soil element or layer which has an initial void ratio e 0 , t 0 is the time corresponding to the completion of primary consolidation, and t is any time beyond primary consolidation. From Equation (8), we can derive the following simplified logarithmic model: where S(t) is the dynamic ground deformation at time t calculated with respect to the reference at the starting time t 0 = 0 (i.e., S(t 0 = 0) ≡0), and a represents the shape parameters of the logarithmic curve. From the literature [26], another model is derived, characterized by an exponential decay of the ground deformations as: (10) where S(t) is the dynamic ground deformation at time t with respect to the reference at the starting time t 0 = 0 (i.e., S(t 0 = 0) ≡0), S max denotes the maximum possible deformation, and S max and b are the parameters of the used model. For comparison, we estimated the ground-deformation time-series from June 2007 to April 2010 at the four points labeled P1, P2, P3, and P4 shown in Figure 13. Then, we obtained the three related best-fit models, which are shown in Figure 15. Then, the quality of the three reconstructed nonlinear curve fittings was evaluated by computing the RMSE of the difference between the InSAR-based (vertical) ground deformation time-series and the obtained best-fit model. We note that the RMSE values between the SBAS-InSAR deformation time-series and the best-fit geotechnical model are smaller than those of the other two models at the four selected points. These results indicate that the geotechnical model outperforms the other two models, also considering the test site in the area of the Shenzhen city.
Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing properties of the reclaimed materials. Finally, we estimated the deformation time-series using the geotechnical model from 2019 to 2020, at all the common pixels between the masked up-down ENV deformation (from 2007 to 2010) map and the vertical S1A deformation (from 2019 to 2020) map. Then, we calculated the grounddeformation velocity by the estimated deformation time-series from 2019 to 2020. To investigate the quality of the ground deformation predictions using the geotechnical model in the whole region, we computed the difference between the vertical S1A deformation velocity and the estimated deformation velocity. The results are shown in Figure 16. This figure evidences that the deformation velocity predicted by the geotechnical model is consistent with and the actual ground deformation velocity recorded in the last year after ten years by the data used for the prediction, obtained by the S1A SAR images. We note that the differences between the deformation velocity predicted by the geotechnical model and the deformation velocity obtained by the S1A SAR data are on average less than 1 mm/year. Accordingly, the results suggest that there is a good agreement between the two deformation velocities. We have also noticed that the deformations are more obvious (lower than −5 mm/year) from 2007 to 2010, while the deformations become very slow from 2019 to 2010 (−2 ~ 1 mm/year), especially in the area highlighted by the red polygons (in Figure 16). Thus, the geotechnical model is reasonable and useful for the deformation prediction of the reclamation area of Shenzhen. Our work, in synergy with the findings of other scholars [3,28,36] in the analyzed area of the city of Shenzhen, is useful to evaluate the risk to the population related to the mutual effects of terrain displacement and present-day sea-level rise (SLR) in coastal regions of highly urbanized cities.  . The plots of the ENV (vertical) deformation time-series relative to the four pixels, labeled as P1, P2, P3, and P4, are shown. The ENV deformation time-series are plotted in pink circles, whereas the corresponding best-fit geotechnical models, logarithmic models, and exponential decay models are plotted by continuous red, green, and blue lines, respectively. The black vertical dotted lines correspond to the present stage (July 2020).
We also estimated the deformation time-series from March 2010 to March 2024 based on the best-fit model parameters of three models at the points P1, P2, P3, and P4, as shown in Figure 15. Compared with the deformation from June 2007 to April 2010, the deformation predicted by the three models from June 2019 to June 2020 is significantly slowed down, and they are generally consistent with the S1A ground deformation measurements. At this stage, we used the independent set of S1A SAR data, collected from June 2019 to June 2020, to check if the measured ground deformation rate in 2020 is consistent with that foreseen by the best-fit models ten years later, based on three years of data collected by the ENVISAT ASAR sensor from 2007 to 2010. As evident, in this case, with a gap of almost ten years, the method proposed in [23] could become unstable and have only a mathematical meaning, so we simply compared the average predicted and measured ground deformation velocities from June 2019 to June 2020. We listed in Table 2 the ground deformation velocities obtained by the ENVISAT ASAR data from 2007 to 2010, those obtained by the S1A sensor from 2019 to 2020, and the ground deformation velocities estimated by the three models from 2019 to 2020. We note that the deformation velocities of S1A and three models are significantly slowed down with respect to the ENV deformation velocities. In terms of the difference between the velocities obtained by three models and S1A deformation velocities, it is evident that the geotechnical model and the logarithmic model are better than the exponential decay model. Furthermore, the differences between the geotechnical model, the logarithmic model, and S1A deformation velocities are in the range of +/−2 mm/year. Moreover, the velocities obtained by the logarithmic model seem to be closer to S1A deformation velocities, especially at points P2, P3, and P4 into the reclamation area. These results indicate that the geotechnical model adopted for the Shanghai city is also suitable in the Shenzhen areas, even though the reclaimed materials are different. Of course, the provided analyses are not conclusive and thus require further efforts and the support of experts of the specific reclaimed-lands project in Shenzhen to relate the model parameters to physical properties of the reclaimed materials.
Finally, we estimated the deformation time-series using the geotechnical model from 2019 to 2020, at all the common pixels between the masked up-down ENV deformation (from 2007 to 2010) map and the vertical S1A deformation (from 2019 to 2020) map. Then, we calculated the ground-deformation velocity by the estimated deformation time-series from 2019 to 2020. To investigate the quality of the ground deformation predictions using the geotechnical model in the whole region, we computed the difference between the vertical S1A deformation velocity and the estimated deformation velocity. The results are shown in Figure 16. This figure evidences that the deformation velocity predicted by the geotechnical model is consistent with and the actual ground deformation velocity recorded in the last year after ten years by the data used for the prediction, obtained by the S1A SAR images. We note that the differences between the deformation velocity predicted by the geotechnical model and the deformation velocity obtained by the S1A SAR data are on average less than 1 mm/year. Accordingly, the results suggest that there is a good agreement between the two deformation velocities. We have also noticed that the deformations are more obvious (lower than −5 mm/year) from 2007 to 2010, while the deformations become very slow from 2019 to 2010 (−2~1 mm/year), especially in the area highlighted by the red polygons (in Figure 16). Thus, the geotechnical model is reasonable and useful for the deformation prediction of the reclamation area of Shenzhen. Our work, in synergy with the findings of other scholars [3,28,36] in the analyzed area of the city of Shenzhen, is useful to evaluate the risk to the population related to the mutual effects of terrain displacement and present-day sea-level rise (SLR) in coastal regions of highly urbanized cities.
Remote Sens. 2020, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/remotesensing deformation velocities estimated by three models from 2019 to 2020 at four pixels labeled as P1, P2, P3, and P4.  Figure 16. Shenzhen area results. The difference between the estimated vertical ground deformation obtained by using the geotechnical model and the S1A vertical deformation products, from June 2019 to June 2020. Figure 16. Shenzhen area results. The difference between the estimated vertical ground deformation obtained by using the geotechnical model and the S1A vertical deformation products, from June 2019 to June 2020.

Conclusions
In this work, we performed an extended investigation to characterize the residual ground displacement of the terrain due to secondary compression stages of unconsolidated ocean-reclaimed lands. In particular, our investigation enhances and extends the outcomes of other previous investigations to another area of interest, namely Shenzhen city, by studying the accuracy of the adopted ground deformation models in a more general context. A comparative analysis involving different sets of independent SAR images has been performed, and the potential and limits of different ground deformation models have been highlighted. We have shed light on the usefulness and applicability of the adopted models to forecast the future evolution of the terrain ground displacement. The results indicate that the geotechnical model used to characterize the deformation phenomena in the Shanghai area can also be successfully used for describing the deformation of other reclaimed areas even though the sub-soil characteristics are not entirely homogenous in the two investigated areas of Shanghai and Shenzhen. The current state of the ground deformations in the area of Shenzhen was obtained using a set of Sentinel-1A SAR data acquired in the last year. Then, the current ground deformation rates were compared with those predicted using sets of ENVISAT ASAR images collected more than ten years ago (from 2007 to 2010) over the same area and already used in investigations performed by other scholars. The results show that the geotechnical model and the logarithmic model are more reasonable in the prediction of the deformation in Shenzhen. Taking the results of Shanghai into account, they allow us to say that the geotechnical model is widely applicable to the description of the deformation in the reclamation area. Globally, the innovation potential of our investigation with respect to previous work [2,19,22,23,[55][56][57] mainly concerns: (i) the corroboration with independent sets of SAR data of the congruence between the measured and the foreseen ground deformation of ocean-reclaimed platforms [2,23]; (ii) the comparative study of different time-dependent models of ground deformation [24,26,48]; (iii) the in-depth analysis of present-day ground displacement in the area of Shenzhen, which has highlighted the presence of newly spotted areas with significant displacement that can evolve over the time and represent a risk to the population [58]. The interrelations existing between the terrain deformation of the Shenzhen area and the present-day sea-level rise (SLR) phenomenon must be seriously taken into consideration. Indeed, even if the rate of the ground deformation of the coastal area is tiny, the results of our investigation reveal that more than a few millimeters of cumulative deformation is expected in the next years and, globally, the process of consolidation of ocean-reclaimed lands continue over an elapsed time of decades. This investigation, entirely based on interferometric synthetic aperture radar data and products, makes it evident that local authorities and scientists have to work more for the achievement of global sustainable development goals.