## 5 Discussion
The map of residuals obtained using the third-order polynomial, although it already gives good results in terms of RMSE and adjusted R\({}^{2}\) (Table 2), shows some relevant boundary effects in the peripheral sectors of the caldera (Figure 10), due to some oscillation in the tail of the curve of the third-order polynomial function (Figure 6b). Such boundary effects are no longer present when using the fourth-order polynomial (Figures 6c and 12). The results obtained with fifth- and sixth-order polynomials (Figure 6d,e, Figures 13 and 14) are comparable to those of the fourth-order polynomial, so there is no need to use these higher orders. Figures 11-14 show that the results of the deformation residuals (second-order anomalies) are very stable in the central sector of the caldera and do not depend on the polynomial function assumed. Figure 18: Plot of normalized vertical ground deformation values (N) vs. distance to the caldera deformation center (R). Comparison between the curve by [PERSON] and [PERSON] [102], which was derived by levelling measurements between 1905 and 1992, 2000–2019 GPS data, and 2016–2021 GPS data [106; 107] used for cross-validation, and our curve, which was derived from Sentinel-1 ortho-vertical DInSAR 2016-2021 data [109]. The different polynomial orders applied to the dataset highlight the complexity of the deformation field and the necessity of balancing between model accuracy and physical interpretation. The fourth-degree polynomial was selected as the best compromise between model complexity and fitting accuracy, since the RMSE and adjusted R\({}^{2}\) values become stable from this order onwards. Higher-degree polynomials only marginally improve model performance but may introduce unnecessary overfitting effects. The deviations from the best fitting curve are significantly non-random (Figure 12). The most prominent lineament is the sharp change between negative and positive anomalies, separated by a roughly NE-SW line slightly east of the Rione Terra and crossing Solfatara crater (Figure 15). This line closely resembles the tectonic lineation connecting La Pietra to the western Agnano plain-bounding faults [26; 110; 111], and most of the seismicity, particularly the largest magnitude earthquakes, appears located along or very close to this line. The clustering of seismicity along the separation line between negative and positive anomalies strongly suggests this is a major structural feature, as does the observed deviation of the deformation pattern from a purely radial one, which is somewhat related to co-seismic displacements. We recall that very shallow seismicity, typical of Campi Flegrei earthquakes, can produce measurable co-seismic displacements from some mm to some cm, even with magnitudes as low as 4. The example of the 2017 Casamicicola (Ischia Island) M = 4.0 earthquake is worth noting; it caused more than 4 cm of subsidence in the hanging wall. The depth of the Casamicicola earthquake was about 2 km, a depth very similar to the Campi Flegrei earthquakes [18; 55; 115; 116]. The detected alignment A in Figure 15, which sharply separates positive and negative vertical deformation anomalies, is then likely to correspond to a main seismic fault zone, which can cause the largest-magnitude earthquakes in the area. It is interesting to note that earthquakes that occurred before 2020 (Figure 7), when the total uplift was lower (65 cm since 2005; Figure 2b) than today, were much more concentrated around this structure. This agrees with the hypothesis that this could be the main seismogenic structure in the area and that other structural features located more distant from the center of uplift have been progressively activated by a progressive enlargement of the critical stress threshold, which is associated with an increase in ground uplift [117; 118; 119]). Such a main structural feature in the area possibly also represents the most likely zone for possible vent opening in a future eruptive scenario (phreatic and/or magmatic).