Before performing micropaleontological analyses, samples were washed and sieved (63 μm) on the sedimentary fraction > 150 μm. The recent benthic foraminifers of the East-Corsica basin were studied in 45 surface samples from the interface cores collected at depths ranging from 7 to 868 m, and 101 benthic foraminifer taxa were identified (Angue Minto’o et al., 2013). The identification of 84 taxa of fossil benthic foraminifer was possible via the analysis of 291 samples from the GDEC-4 borehole drilled at a water depth of 492 m in the East-Corsica margin covering the last 550,000 years (Angue Minto’o et al., 2016).

Oxygen stable isotope measurements were performed on specimens of planktonic foraminifera species Globigerina bulloides and G. ruber (white) from the 250–315 μm size fraction, Neogloboquadrina pachyderma (dextral) from the 200–250 μm size fraction, and on the epifaunal benthic foraminifera Cibicides wuellerstorfi, Cibicidoides pachyderma and Cibicidoides kullenbergi found in the > 150 μm size fraction (Toucanne et al., 2015).

In this study, MatLab generic functions are used for the computations. We call M₀(i,j), the initial matrix of species to be consider in the analysis. M₀(i,j) consists of relative abundances of the all benthic foraminifers (101 species) identified in the surface samples, where i is the site and j the species. With the aim of obtaining qualitative results, the reduction of the number of species of M₀(i,j) is made by eliminating species with a median equal to 0. Because Principal Component Analysis (PCA) is based on correlation analysis and hence variation quantifications, consequently, on 101 recent benthic foraminifer taxa identified, only 33 taxa were retained and are listed in Table 1. For these 33 species, the relative abundance, at one site, is always calculated using the total number of individuals per site based on the 101 species. This allowed maintaining independence between the frequencies retained, i.e. the sum of the species frequencies per site does not equal 1, and the dependent and untreated variable is “the other species”

PCA is based on the covariance matrix (M_c) estimated on the basis of the species abundance matrix by site:

The eigenvalues and eigenvectors of the matrix M_c are computed to establish, respectively, the weight (variance) associated with each component and the coefficients of the principal components. Thus, the matrix of the eigenvectors C_jk, and the coefficients of the principal components are calculated as well as the diagonal matrix of the eigenvalues D_jk. The new matrix C_jk verifies the following equality:

The principal components are obtained by a “rotation” of the assembly matrix (M_ij) along the principal vectors described by the matrix of eigenvectors or coefficients of the components

It is important to keep in mind that the matrix Z_ik contains the same amount of information as the initial abundance matrix M_ij from which Z_ik has been derived. The information is organised in a new way, the components being independent of one another. Hence, at each site i, Z_ik provides the local value of component k, while M_ij provides the local abundance of species j, the number of species and of components being the same (33).

The matrix C_jk co-relates the principal components and the species. Even if it is computed from the data, it contains an estimation of a relation that is supposed to be valid for any assemblage (not only the ones in the data). C_jk is the same for all sites (observed or noted) and is the base of the generalization.

At each new site,

We can evaluate the possibility of obtaining a model linking the species assemblages and the depth of the sites.

For each component Z_k, one can define a depth model H_k obtained by linear regression for the depth H:

For models consisting of several components, the method applied is the same. A depth model is obtained by linear regression using the least squares method. Here the most representative component of the variance is used. Thus, a depth model H_n can be defined, where n is the number of components used (from the highest to the lowest):

The performance of the model H_n can be associated with the A_n function,

The value A₆ will give the correlation level of a regression model based on the six most important components in terms of assemblage variability.

A freedom factor (λ_n) of a model consisting of principal components n with respect to the points number (P) of calibration/validation is defined by the following formula:

A freedom factor is defined to quantify the robustness of the models based on the number of its number of freedom degrees. For a model based on n principal components, its freedom factor (L_n) is quantified by the number of data points used for calibration (P) and the number of degrees of freedom of the used model (n + 1, if n PCs are used).

If this factor is close to the unit, the model is robust in terms of degrees of freedom. If the factor is close to zero, the model is not robust, and there are too many degrees of freedom compared to the size of the data available and the basis of its construction. In our study, P = 44.

An performance indicator of the model could be defined as the “good” compromise between the correlation level and its performance or factor of freedom. The quality of a model H_n described by the function Q_n is defined as follows:

After a ranking based on their depths, the 44 sites are divided into two groups. The first one is composed of one site out of two and it is called “calibration”. This group includes variables with an index c. The second group consists of the other sites and it is called “validation”. This validation group is characterized by variables with an index v. The coefficients a_k are calculated by regression on the basis of the following system of equations:

The modelled values will then be:

The calibration coefficients are then used for the validation data in order to obtain the simulated data H_m,v:

The results of the PCA show that a very important part of the variance (> 90%) is contained in the first five eigenvalues, or principal components (Fig. 1). The contribution of each species to the two main components can be visualized by the respective contributions of each species (Fig. 2). It may be noted in Fig. 2 that the species that contribute the most to the expression of the first component are M. barleeanus (21) and U. mediterranea (33). This component is modulated (attenuated) by the presence of C. carinata (11), R. globularis (26), and Q. duthiersi (23). The second component, however, is expressed by the “competition” between C. carinata and R. globularis (Fig. 2).

The correlation levels between the depth of the sites and the principal components associated with the assemblages are presented in Fig. 3. The function B is represented by the blue bars that indicate the correlation level of each component with the depth. The function A is represented by the red curve. This latter indicates the correlation level reached by regression between a set of N components and the depth. Thus, the component 33, which is the first in importance, is the one that shows the best correlation with the depth. This component alone allows having a correlation level higher than 95% (Fig. 3). The other components are individually weaker. However, they allow increasing the total correlation to more than 99% if all 33 components are used to reproduce the 44 depths.

The best model is not necessarily the one that allows obtaining the best correlation between true values and model values estimated on these same true values. The number of degrees of freedom associated with the calibrated models and the number of validation points must be taken into account, because the models are calibrated on the validation points. Thus, all models H_k have two degrees of freedom: the coefficients a and b. The models H_n have n + 1 degrees of freedom. By absurdity, a model whose number of degrees of freedom is equal to the number of calibration/validation points will always present a correlation of 100%.

The representation of the function Q_n (Fig. 4) shows that the use of a small number of components can be advantageous in terms of robustness of the chosen model.

The evaluation is made for models built with a number of variable components (between 1 and 6). The model with all the data (blue curve) serves as a reference (Fig. 5). Calibration 1 (cal 1) and validation 1 (val 1) are obtained as described in the paragraph related to calibration and validation of models: the first set is for calibration and the second set for validation. The models cal. 2 and val. 2 are obtained by reversing the use of the sets (Fig. 5).

For cal. and val. 1, we obtain a calibration model that is slightly less efficient than the reference one, while the validation is generally better (Fig. 5). This means that the validation set is closer, in terms of assemblages, to what can be explained by the depth. Overall, the correlation improves with the number of components used. This improvement is increasing for the reference and for the calibration (Fig. 5). For the validation model, the two minor contribution components (2 and 5) reduce the performance of the model. When we reverse the use of sets, we fall back on a reverse result, and more “classic”. Calibration, done on less data than the reference, is more efficient (but with a lower freedom factor); validation is less well than calibration and reference.

The choice seems to be made between 1, 4, or 6 components, even if it remains conceivable to use only the useful components. Correlations are performance indicators. It is therefore important to visualize the correlations presented in Figs. 6 and 7, in which a loss of linearity can be observed from 700 m depth with the error margins of 43 m and thus 86 m in 2σ.

Principal component analysis allows us to “focus” the variance of the assemblages on a very small number of variables (components), which is very useful in order to reduce the number of freedom degrees of the models obtained by regression.

Out of 33 components (as much as the selected species), we will keep four components (28, 30, 31 and 33), which means that the regression model is based on the estimation of five parameters.

The model obtained made it possible to have a 97.13% correlation on all the sites. The calibrated model on 22 sites (performance of 96.85% in self-application) displays a correlation of 97.72% on the validation set. The performance of the chosen model is illustrated in the Fig. 8.

The structure of the model is presented in Fig. 9. It can be noted that the final contribution to explaining the depth variability is very distinct from the distribution of the coefficients associated with each species. This is because the relative abundance of each species has not been standardized. Thus, there are two species that contribute significantly: M. barleeanus and U. mediterranea. M. barleeanus explains the variations of depth at shallow depths and U. mediterranea explains the depth beyond 100 m.

4.6.1 Numerical model

The application of the model H(i) = 109.6103 + C₁ × F1(i) + C₂ × F₂(i) + … + C₃₃ × F₃₃(i) (where C_i corresponds to the coefficient of species i and F_i to its frequency) on fossil benthic foraminifers species allowed us to obtain the curve of paleo-depths variation. The reconstructed paleo-depths range from 622 ± 86 m to 21 ± 86 m. In general, the variation of these paleo-depths correlates quite well with the variations of the sea level (Fig. 10). A global trend marked by an increase in depth during periods of high sea level is observed. Interglacials are characterized by the removal (or melting) of glaciers that conduct to an increase in sea level (Dorale et al., 2010) and therefore to an increase in the water level above the seabed. This could justify this rise of the paleo-depths during these warm climatic periods. However, these paleo-depth variations are characterized by very large amplitudes.

Based on the fact that these modelled amplitudes are affected by (1) local effects related to the morphology of the basin, which is shallow and semi-marine (close to the sources), (2) changes in trophic conditions at the bottom strongly influencing the variations of benthic microfauna assemblages, (3) uncertainties related to the calculation method and the bathymetry. Normalization from 0 to 120 m was thus made on the values of paleo-depths and the sea level variation curve obtained was compared with the eustatic curve established by Rohling et al. (2009) and smoothed in the same way (Fig. 11). This results in a fairly good correlation between these two curves.