M A Pais

A L Morozova

N Schaeffer

Variability modes in core flows inverted from geomagnetic field models

The flow of liquid metal inside the Earth's core produces the geomagnetic field and its time variations. Understanding the variability of those deep currents is crucial to improve the forecast of geomagnetic field variations, which affect human spacial and aeronautic activities. Moreover, it may provide relevant information on the core dynamics. The main goal of this study is to extract and characterize the leading variability modes of core flows over centennial periods, and to assess their statistical robustness. To this end, we use flows that we invert from two geomagnetic field models (gufm1 and COV-OBS), and apply Principal Component Analysis and Singular Value Decomposition of coupled fields. The quasi geostrophic (QG) flows inverted from both geomagnetic field models show similar features. However, COV-OBS has a less energetic mean and larger time variability. The statistical significance of flow components is tested from analyses performed on subareas of the whole domain. Bootstrapping methods are also used to extract significant flow features required by both gufm1 and COV-OBS.

Three main empirical circulation modes emerge, simultaneously constrained by both geomagnetic field models and expected to be robust against the particular a priori used to build them (large scale QG dynamics). Mode 1 exhibits three large vortices at medium/high latitudes, with opposite circulation under the Atlantic and the Pacific hemispheres. Mode 2 interestingly accounts for most of the variations of the Earth's core angular momentum. In this mode, the regions close to the tangent cylinder and to the equator are correlated, and oscillate with a period between 80 and 90 years. Each of these two modes is energetic enough to alter the mean flow, sometimes reinforcing the eccentric gyre, and other times breaking it up into smaller circulations. The three main circulation modes added together to the mean flow account for about 70% of the flows variability, 90% of the root mean square total velocities, and 95% of the secular variation induced by the total flows.

Direct physical interpretation of the computed modes is not straightforward. Nonetheless, similarities found between the two first modes and time/spatial features identified in different studies of core dynamics, suggest that our approach can help to pinpoint the relevant physical processes inside the core on centennial timescales. 1

Introduction

Time variation of the main geomagnetic field on timescales smaller than about a couple of centuries (secular variation, SV) have been scrutinized for deeper understanding of the Earth's core dynamics. New insights have been gained both from numerical simulation of the Navier-Stokes, energy and induction equations inside the core [see e.g. [START_REF] Jones | Planetary Magnetic Fields and Fluid Dynamos[END_REF] and inversion of magnetic field models fitting observatory and satellite observations [e.g. [START_REF] Holme | Large-scale flow in the core[END_REF]. Although the parameter regime characterizing the Earth's short timescale dynamics is not yet attainable in three-dimensional numerical simulations, different schemes have been used to extrapolate results. These include derivation of scaling laws and/or the use of simplifying approximations that allow the use of more Earth-like parameters.

The Quasi-Geostrophic (QG) approximation is such an example, since it reduces the system dimensionality from three to two (3D to 2D) making it possible to investigate a range of Ekman (ratio of viscous to Coriolis forces) and Lundquist (ratio of Alfvén wave to diffusion timescales) numbers that approach the Earth's core conditions. Convective quasi-geostrophic rolls are known since the work by [START_REF] Roberts | On the thermal instability of a self-gravitating fluid sphere containing heat sources[END_REF] and [START_REF] Busse | Thermal instabilities in rapidly rotating systems[END_REF] and recent studies show that for timescales characteristic of SV and large scale flows, elongated structures along the z-rotation axis will not be destroyed by Lorentz (magnetic) forces [Jault, 2008[START_REF] Gillet | Rationale and geophysical evidence for quasi-geostrophic rapid dynamics within the earth's outer core[END_REF]. One must however keep in mind the limitations of the QG approximation, mainly that it does not allow for anti-symmetric flows (the ones that cross the equator) which could in reality be produced by any non-symmetric forcing.

From the point of view of flows computed from inversion of geomagnetic field models, the QG approximation is one example of constraint used to mitigate the intrinsic non-uniqueness characterizing the inverse problem. Other constraints relying on different dynamical assumptions or empirical reasons have also been used (e.g. purely toroidal flow, steady flow, tangential geostrophy, steady flow in a drifting reference frame, helical flow, and tangential magnetostrophy). An account on these previously used constraints can be found in [START_REF] Holme | Large-scale flow in the core[END_REF] and [START_REF] Finlay | Short Timescale Core Dynamics: Theory and Observations[END_REF]. In the following, we will restrict to QG flows due to their importance in SV dynamics [START_REF] Gillet | Rationale and geophysical evidence for quasi-geostrophic rapid dynamics within the earth's outer core[END_REF]. In addition, although the flow is inverted on the core surface, it has a unique prolongation inside the core. QG flows inverted from different geomagnetic field models have been previously computed that explain the observed SV [Pais andJault, 2008, Gillet et al., 2009]. In these studies, a large scale spatial structure was identified and described as a large eccentric anticyclonic gyre. It interestingly gathers into one single structure the well-known westward drift under the low latitude Atlantic region, and the high latitude jet under the Bering Sea and close to the tangent cylinder (TC, cylinder coaxial with the Earth's rotation axis and tangent to the inner core surface) [see e.g. Figure 3 in [START_REF] Holme | Large-scale flow in the core[END_REF]. This structure seems to be more prominent in inversions derived from recent satellite geomagnetic field models, suggesting some time variability on relatively short timescales characteristic of SV. Other structures are also retrieved, a large scale cyclonic vortex beneath the Pacific Hemisphere and smaller scale vortices dominantly cyclonic under the Pacific Hemisphere and anticyclonic under the Atlantic Hemisphere. [START_REF] Gillet | Ensemble inversion of time-dependent core flow models[END_REF] pointed out that the smaller scale vortices may not in fact be well resolved, since they strongly depend on the induction effects that involve the small scale main field (above harmonic degree 13-14), which are difficult to infer reliably from surface and satellite observations. In this respect, [START_REF] Finlay | Core surface magnetic field evolution 2000-2010[END_REF] stress the fact that the spectral slope above degree 12 depends significantly on the modeling regularizations in space. But there has not yet been a dedicated study to establish the significance of the flow structures. In this paper, we study the time variability of the inverted flow, and we assess their significance with statistical tools largely used in climate studies and oceanography.

Recently, a new geomagnetic field model has been computed, COV-OBS, which covers the period of historical geomagnetic field data provided by the international network of observatories , but also includes the recent period of satellite missions with a dense coverage of high quality data [START_REF] Gillet | Stochastic modeling of the earth's magnetic field: Inversion for covariances over the observatory era[END_REF]. It avoids employing the regularized least squares inversion approach used to compute the gufm1 model [START_REF] Jackson | Four centuries of geomagnetic secular variation from historical records[END_REF], where roughness in the spatial and in the temporal domains is minimized when fitting the data to a truncated spherical harmonic expansion and using low-order splines in the time domain. Claiming that such regularizations are mainly ad hoc and do completely constrain the model above some harmonic degree (smaller at older epochs, larger at recent ones), [START_REF] Gillet | Stochastic modeling of the earth's magnetic field: Inversion for covariances over the observatory era[END_REF] preferred to use a priori information on the statistical properties of the geomagnetic field, which they consider to be relatively well known, and follow a stochastic modeling approach (although their approach is not completely free from some imposed temporal constraints, since the COV-OBS model is expanded on a spline basis with a 2-year knot spacing). The Gauss coefficients of the geomagnetic field are treated as resulting from a stationary process, their statistical properties being condensed into a large prior covariance matrix with variances and time correlations for all coefficients. Although the two models gufm1 and COV-OBS fit the shared data at the same level for the common period , the different intrinsic characteristics are expected to distinguish flows inverted from one and the other.

In this study, we use QG flow models inverted from both gufm1 and from COV-OBS. The computation of these flows, which is outlined in section 2, is explained and discussed in the appendix A. Our main purpose is to test the possibility to describe the global circulation of the core flow in terms of a small number of patterns of variability, much like the climate variability can be described as a combination of patterns like the ENSO (El Niño-Southern Oscillation), PNA (Pacific -North America teleconnection pattern), NAO (North Atlantic Oscillation), etc [e.g. [START_REF] Barry | Synoptic and Dynamic Climatology[END_REF]. We apply Principal Component Analysis (PCA) tools, of standard use in the geophysical sciences of meteorology and oceanography, where we find the main references [e. g Preisendorfer, 1988, Von Storch and[START_REF] Storch | Analyses of Climate Variability-Applications of Statistical Techniques[END_REF]. These are nonparametric methods that do not assume any particular statistical distribution for the analyzed data or any particular dynamical equation for the flow. In geomagnetic studies, PCA has been previously applied to space-time data grids in studies of geoelectromagnetic induction [START_REF] Fujii | The 3D electromagnetic response of the Earth to ring current and auroral oval excitation[END_REF]Schultz, 2002, Balasis and[START_REF] Balasis | Empirical orthogonal function analysis of magnetic observatory data: Further evidence for non-axisymmetric magnetospheric sources for satellite induction studies[END_REF], and using ionospheric data [START_REF] Matsuo | Modes of high-latitude electric field variability derived from DE-2 measurements: Empirical Orthogonal Function (EOF) analysis[END_REF][START_REF] Zhang | Ionospheric symmetry caused by geomagnetic declination over North America[END_REF]. The PCA approach has several limitations, as explained by [START_REF] Richman | Rotation of principal components[END_REF]. These may be due to domain shape dependence through the space orthogonality condition imposed, to sensitivity to the extent of spatial and time domains considered and to incomplete separation of consecutive modes. Also, there may be no correspondence between PCA-computed and physical modes. Nevertheless, PCA variability modes generally allow for a complexity reduction and, in certain cases, can provide information on the underlying dynamics [e. [START_REF] North | Empirical Orthogonal Functions and Normal Modes[END_REF]]. They have accordingly been extensively used in data analysis. Besides PCA, singular Value Decomposition (SVD) of coupled fields is also applied in this study to identify pairs of coupled spatial patterns explaining the covariance between gufm1 and COV-OBS [see e.g. [START_REF] Bjornsson | A Manual for EOF and SVD Analyses of Climatic Data[END_REF]. At the same time, we propose to carry out significance tests that can decide which are the modes carrying relevant information on the system variability. The obtained results receive statistical support from analyses performed on subareas of the whole domain. Bootstrapping resampling methods are used to assess the statistical robustness of particular flow features. In the end, a physical interpretation of the most important/significant modes is sought.

The flow data to be analyzed is described in section 2. The basis of the PCA method and the more technical aspects that are used in this study are introduced in section 3. Results are discussed in section 4. They include results obtained when analyzing flows inverted from gufm1 and from COV-OBS separately, and results obtained by extracting correlated information contained in both of them. Section 5 contains a discussion on the identified structures and their possible dynamical interpretation. Finally, the last section summarizes the main conclusions in this study.

Inverted QG flow models

Two core surface flow models were computed, using regularized inversion, from the timevarying geomagnetic field models for the observatory era gufm1 (1840 -1990) [START_REF] Jackson | Four centuries of geomagnetic secular variation from historical records[END_REF] and COV-OBS (1840-2010) [START_REF] Gillet | Stochastic modeling of the earth's magnetic field: Inversion for covariances over the observatory era[END_REF]. These two geomagnetic field models employ cubic B-splines to model the time variation, with 2.5-year and 2-year knot spacing, respectively (but data are supplied twice the knot spacing). Our flow inversions generate a flow snapshot for each year in the geomagnetic field model period, and the whole set of flows is denoted flow guf m1 and flow COV -OBS , respectively. The whole inversion procedure is explained in detail in the appendix A.

In the inversion, the flows were assumed to be columnar in the whole volume and incompressible. This implies that they can be completely retrieved from a single scalar pseudo streamfunction ξ, symmetric with respect to the equatorial plane, which in turn can be inverted from geomagnetic field models. At the core surface, the relation between the core surface flow u and ξ at a certain epoch is [START_REF] Schaeffer | Quasigeostrophic model of the instabilities of the Stewartson layer in flat and depth-varying containers[END_REF][START_REF] Pais | Quasi-geostrophic flows responsible for the secular variation of the earth's magnetic field[END_REF][START_REF] Amit | Differences between tangential geostrophy and columnar flow[END_REF]:

u = 1 cos θ ∇ H ∧ ξ(θ, φ) r + sin θ R c cos 2 θ ξ N Z (θ, φ) φ (1)
where (r, θ, φ) are spherical coordinates, R c is the core radius, and ξ N Z is the non-zonal (φdependent) part of ξ. The first term in (1) has exactly the same expression as tangential geostrophic (TG) flows computed at the CMB [Chulliat and Hulot, 2000, e.g.], where ξ would represent the geostrophic pressure to within a constant factor, p geo = -2ρΩξ (with ρ the core fluid density and Ω the Earth's rotation rate). For columnar flows touching the CMB, if they conserve mass they no longer follow isocontours of ξ because of the second term in (1). This term does not contribute to the zonal flow, is everywhere latitudinal and dependent on co-latitude. It is required for fluid mass conservation, leading to equal contributions to upwellings and downwelling from the two terms [see e.g. [START_REF] Amit | Differences between tangential geostrophy and columnar flow[END_REF]. As also discussed in [START_REF] Amit | Differences between tangential geostrophy and columnar flow[END_REF], the ξ field is uniquely determined from inversion of geomagnetic field models. Close to the equator, ξ N Z (θ, φ) ∼ cos 2 θ in order for the flow to remain finite. [START_REF] Canet | Hydromagnetic quasi-geostrophic modes in rapidly rotating planetary cores[END_REF] use other pseudo streamfunction Ψ(θ, φ) to study hydromagnetic quasi-geostrophic modes in planetary interiors, the two scalar functions being related through

Ψ N Z = -R c cos θ ξ N Z and ∂ θ Ψ Z = -R c cos θ ∂ θ ξ Z (
where the Z superscript denotes the zonal component).

The assumption of equatorially symmetric flow due to columnar convection does not hold inside the tangent cylinder. Different regimes can be envisioned there, depending on the dominant force balance believed to exist in that region. In particular, this can lead either to independent columnar convection in the North and South hemispheres [START_REF] Pais | Quasi-geostrophic flows responsible for the secular variation of the earth's magnetic field[END_REF] or thermal wind convection [e.g [START_REF] Aurnou | Experiments on convection in Earth's core tangent cylinder[END_REF]] that can change the sign of circulation close to the inner core and the outer core boundaries.

Because of this uncertainty, and for convenience, we use the same assumptions inside and outside the TC for computing u. Hence, to guarantee that flows computed for this small region (that usually tend to have high root mean square velocity u rms ) will not affect the PCA modes, the 20 • caps around the North and South CMB poles were excluded from the analysis as explained in the following section.

Principal Component Analysis of QG flows

We follow the implementation of the method as explained in [START_REF] Bjornsson | A Manual for EOF and SVD Analyses of Climatic Data[END_REF], [START_REF] Hannachi | Empirical orthogonal functions and related techniques in atmospheric science: A review[END_REF]. The field ξ, closely related to a geostrophic pressure especially at high latitudes, can be used to characterize the spatial flow features and its time variability over some time-period, with a similar role to sea-level pressure in studies of variability of atmospheric currents [e.g. [START_REF] Hannachi | Empirical orthogonal functions and related techniques in atmospheric science: A review[END_REF]. In short, we will analyze the ξ field into a sum of orthogonal functions in space, multiplied by coefficients which are uncorrelated functions of time, by applying a linear transformation to the dataset that concentrates as much of the variance as possible into a small number of terms in the expansion.

Empirical Orthogonal Functions and Principal Components

The spatio-temporal information retrieved from inverted flow models is gathered into matrices:

X =      X 11 X 12 • • • X 1Np X 21 X 22 • • • X 2Np . . . . . . . . . . . . X Ne1 X Ne2 • • • X NeNp      (2)
with

X ij ≡ X(t i , r j ) = √ sin θ ξ(t i , θ j , φ j ) - 1 N e Ne i=1 ξ(t i , θ j , φ j ) . (3)
In a few words, X ij is the (scalar) value of the ξ pseudo-streamfunction at the core surface point r j = r(R c , θ j , φ j ) at time t i , once the time-average has been subtracted. It is also weighted by a latitude-dependent factor √ sin θ to compensate for the fact that in a regular grid on a spherical surface the number of points to cover a given area increases with latitude. The index i takes values from 1 to the total number of epochs, N e , and the index j takes values from 1 to the total number of grid points, N p . For flow guf m1 and flow COV -OBS , a set of 5 • × 5 • latitude/longitude grids of the scalar function ξ were computed, one for each year, covering the CMB region under study defined by 20 • < θ < 160 • for all values of longitude. We exclude the polar caps because, as referred above, the computed QG flows are not reliable inside the tangent cylinder (TC).

The Empirical Orthogonal Function/Principal Component analysis (EOF/PCA) relies on the description of the space-time data in terms of decorrelated modes. They are the eigenvectors of the covariance (or variance-covariance) matrix C X = X T X: each element (C X ) ij = X li X lj = Ne l=1 X(t l , r i )X(t l , r j ) is the discrete form of the (temporal) covariance between the data values at points r i and r j . In search for decorrelated modes, we need to transform C X in a diagonal matrix, for which the standard procedure is to transform the data matrix X into Y = XP (4)

where P is the orthogonal matrix of eigenvectors of C X , each eigenvector p i forming a column.

The covariance matrix of Y is then

C Y = P T C X P = Λ (5)
and is diagonal with non-zero elements λ i , ordered from the highest to the lowest. The transformed data-matrix Y is now expressed in terms of decorrelated modes, as was intended. From (4) the reconstructed data matrix yields

X = YP T , (6)
which can be written

X = y 1 p T 1 + y 2 p T 2 + • • • + y Np p T Np = Np k=1 y k p T k , (7)
i.e., as a linear combination of all eigenvectors of C X , the coefficients of the expansion being the different columns of matrix Y.

Here, and according to [START_REF] Bjornsson | A Manual for EOF and SVD Analyses of Climatic Data[END_REF], the following notation is used:

• each column p k of P is an Empirical Orthogonal Function (EOF) and represents a certain spatial function;

• each column y k of Y is a Principal Component (PC) and represents a time series.

This decomposition presents the further advantage of data reduction, since there is no practical need to keep the whole set of N p EOF s (the same number as grid points) in the analysis. The different eigenvalues λ i of C X , determine which modes are to be kept in the expansion and in practice only the first few ones are important, those that explain the highest percentage of data variance as given by

f i = λ i   Np k=1 λ k   -1 . ( 8
)
In terms of describing the data covariance matrix, upon which the whole framework is based, the EOFs provide a simplified and (hopefully) compressed way to decompose this matrix. From (5) it follows

C X = Np k=1 λ k p k p T k , (9)
with only the first modes required. We note that some studies apply PCA using correlation instead of covariance matrices, by normalizing each time series to unit standard deviation [see e.g. [START_REF] Bretherton | An Intercomparison of Methods for Finding Coupled Patterns in Climate Data[END_REF]. Here, such a choice could enhance small, non-meaningful signal and was not adopted.

Singular Value Decomposition of coupled fields

Singular Value Decomposition (SVD) of coupled fields has been widely used in meteorology and oceanography, to study two combined data fields of different physical quantities such as e.g. Sea Level Pressure (SLP) and Sea Suface Temperature (SST). In this study, it is used as a framework to isolate important coupled modes of variability between the two time series of core flows, flow guf m1 and flow COV -OBS , representing the same physical quantity but derived from geomagnetic field models computed using different a priori on the data and the model itself. We argue that by applying SVD to the correlation between the two flows the most robust features can be recovered, those that do not depend specifically on the a priori used in one or the other geomagnetic field models.

We follow the implementation of the method as explained in [START_REF] Bretherton | An Intercomparison of Methods for Finding Coupled Patterns in Climate Data[END_REF], [START_REF] Bjornsson | A Manual for EOF and SVD Analyses of Climatic Data[END_REF] and [START_REF] Hannachi | Empirical orthogonal functions and related techniques in atmospheric science: A review[END_REF]. Starting from two data matrices X 1 and X 2 for flow guf m1 and flow COV -OBS respectively, built as explained above, the temporal cross-variance matrix is constructed,

C X 1 X 2 = X T 1 X 2 .
As will be made clear in the following, it makes no difference for the spatial patterns, the temporal functions or the percentage of cross-covariance explained, that this covariance matrix or its transpose, C X 2 X 1 , be considered. The SVD general matrix operation is applied to C X 1 X 2 :

C X 1 X 2 = UΛV T , (10)
where U and V are both orthogonal matrices, the former having in columns the eigenvectors of

C X 1 X 2 C T X 1 X 2 and the later having the eigenvectors of C T X 1 X 2 C X 1 X 2 .
As to Λ, it is a diagonal matrix with singular values λ i in its diagonal, ordered from the highest to the lowest. These λ i are the (positive) square roots of the eigenvalues of C T X 1 X 2 C X 1 X 2 , the same as for C X 1 X 2 C T X 1 X 2 . Linear transformation of data matrices X 1 and X 2 with the orthogonal matrices U and V respectively, makes the new cross-correlation matrix diagonal, C Y 1 Y 2 = Λ, where Y 1 = X 1 U and Y 2 = X 2 V. The original data matrices are reconstructed as

X 1 = Y 1 U T X 2 = Y 2 V T , (11)
analogously to (6). It is nonetheless worth noting that U and V integrate information on both flow guf m1 and flow COV -OBS , contrary to the transformation matrix P in section 3.1. As a result, spatial patterns obtained from U for instance are different from EOFs obtained by applying PCA to X 1 . Also, the expansion coefficients time series in (11), namely Y 1 and Y 2 , are different from the PCs referred in (3.1).

In analogy with (9), the column vectors u i and v i of the transformation matrices U and V, respectively, can be used to decompose the matrix C X 1 X 2 :

C X 1 X 2 = Np k=1 λ k u k v T k . ( 12
)
This expression resorts the fact that each pair k of coupled modes describes a fraction of the cross-covariance between the two fields, the fraction being given by f k (see 8).

Homogeneous and heterogeneous correlation maps

Besides direct representation of PCA or SVD modes in the form of maps of the EOFs and time plots of corresponding PCs, correlation maps are a means to highlight in the analyzed spatial region, those regions that are evolving in time in a correlated manner. In each point of the space grid, a correlation coefficient is computed between a PC time function and the 'observed' time-evolving streamfunction in that same point. For homogeneous correlation maps, the concept is very similar to the one leading to different EOFs and correlating the different PCs with the original data gives rise to spatial structures that are the same as the corresponding EOFs. Nonetheless, it involves looking for strong correlation instead of strong covariance and, as such, it is normalised to local variance and the values in each grid point are between -1 and +1. In this way, regions of low amplitude variations that show up as faint regions in the charts of EOFs, will show up as strong as regions of high amplitude variations, if they belong to the same EOF structure. Heterogeneous correlation maps are meaningful in SVD analysis to show how PCs extracted from a certain flow model are correlated with the other flow model. If the identified SVD modes are to be given some physical meaning, it is required that a certain PC correlated with the two flow models gives similar structures.

From studies of climatic patterns analyses, we borrow standard tests used to assess the robustness of variability modes. The relevant quantities are introduced in appendix B.

Results of statistical analyses

The mean flows shown in Fig. 1 are computed by averaging the flow over the corresponding time-period: 1840 -1990 for flow guf m1 and 1840 -2010 for flow COV -OBS . There are no significant differences if we consider the shorter 1840 -1990 period for COV-OBS. The root mean square (rms) velocity for the surface mean flow, u rms , is stronger for flow guf m1 (12 km/yr) than for flow COV -OBS (10 km/yr), basically due to a stronger anticyclone under the Pacific Hemisphere and a stronger radial jet under the Eastern Asian continent, from high to low latitudes. The two different kind of charts in Fig. 1 illustrate the fact that the flow very closely follows contours of the ξ streamfunction at medium-high latitudes. The rule to keep in mind when interpreting ξ-contour charts, from eq. 1, is that the fluid circulates anti-clockwise around centers of positive ξ and clockwise around centers of negative ξ.

The EOF/PCA tools can sort out the time variability associated to the main structures in the mean flow, and other structures which average out during the inspected period. In the following, we will be showing results for the analysis carried out over the computed flows, using tools described in section 3 and appendix B.

PCA applied to gufm1 and COV-OBS, separately

The first five PCA modes account for about 80% of the flow total variability (see Table 1). The EOF patterns of the first three modes are shown in Figure 2. Results for COV-OBS considering the two time periods are similar. Note however that the first 3 modes represent a larger fraction of the signal for the shorter period. Besides, the 4th mode may not be completely separated from the 3rd when considering the longer time period. Comparing COV-OBS and gufm1, the first mode explains very similar variances. Using North's criterion (eq. 22) to detect mode degeneracy, modes 4 and 5 in COV-OBS (1840 -1990) and modes 3 and 4 in COV-OBS (1840 -2010) are not completely separated (see Figure 3), meaning that

Fraction of variance,

f i [%] (δf i [%]) mode-1
mode-2 mode-3 mode-4 mode-5 COV-OBS (1840 -2010) 34.1 (4.1) 16.6 (2.0) 11.0 (1.3) 10.3 (1.2) 7.5 (0.9) COV-OBS (1840 -1990) 34.4 (4.5) 18.3 (2.4) 12.7 (1.7) 8.6 (1.1) 7.7 (1.0) gufm1 (1840 -1990) 33.3 (4.3) 22.9 (3.0) 16.2 (2.1) 7.9 (1.0) 6.6 (0.9)

Table 1: Standard PCA applied to flow guf m1 and flow COV -OBS separately. Percentage of variance explained by each mode (100f i , using eq. 8). In parenthesis, the δf i error due to δλ i as given by eq. 21. they may describe two aspects of a common structure. A simple illustration of this effect is found in the example of a propagating wave that can be decomposed into two spatial patterns space-shifted by fourth of a wavelength, multiplied by two sinusoidal functions time-shifted by fourth of a period. One unique structure (a propagating wave) would then appear in the PC analysis decomposed into two modes with exactly the same f value. In the case of the above-mentioned degenerate modes, they should be considered together.

As a test for subdomain stability, the whole CMB domain was subdivided into two longitudinal hemispheres, and EOF/PC modes were recalculated for each of them. The chosen meridian for the separation goes through 70 • E and the two resulting hemispheres have longitudes 70 • E to 250 • E (Pacific Hemisphere, PH) and -110 • E to 70 • E (Atlantic Hemisphere, AH). Results are shown in Tables 2 and3. The first five modes are still non-degenerate when computed independently for each hemisphere and according to North's criterion (eq. 22) (see Table 2). While the spatial and temporal descriptions of modes 4 and 5 can be different depending on if a global or an hemispherical grid of data values is used, the first three modes are recovered under AH with very close characteristics as in the global grid. The first three variability modes are not so well recovered using only data under the PH, especially for mode 2 (see Table 3). Overall, when applying PCA to gufm1 and COV-OBS separately, some significant differences come out: (i) mode 2 is more important to explain the variability of flow guf m1 than of flow COV -OBS , possibly due to a less clear expression of mode 2 under the PH; (ii) also, differences between recovered modes in PH and AH hemispheres are more important for flow COV -OBS than for flow guf m1 ; (iii) most importantly, corresponding modes in flow guf m1 and flow COV -OBS exhibit significant spatial and temporal differences in general, as shown by relatively low values of congruence and correlation coefficients in square brackets in Table 4 (see also Fig 2).

SVD of coupled gufm1 and COV-OBS

In order to identify structures correlated in both flows, we apply SVD of coupled fields to flow guf m1 and flow COV -OBS for the common period 1840-1990 (see Table 4).

In Figure 4 we represent the contours of streamfunction ξ for the three first coupled modes. As all modes show important spatial structures at medium-high latitudes, where visualization of the flow through ξ contours is quite straightforward, the representation adopted uses these contours projected onto the equatorial plane as seen from the North pole. Each pattern EOF i contributes to the global flow with the flow circulation retrieved directly from Fig. 4 in epochs of positive PC i and opposite circulation in epochs of negative PC i . We easily recognize common main structures in the pair of EOFs for the first two modes, namely three large rolls at medium-high latitudes for mode 1, two anti-cyclonic under the AH and one cyclonic under the PH, and an important zonal component around the TC in mode 2. Mode 2 also presents a strong zonal component at the equator, which is not apparent in Figure 4 because of the chosen projection (see however Fig. 9 and the Discussion in section 5). Special care must be taken in interpreting the EOF 2 pattern in Figure 4. Indeed, the streamfunction ξ has an important negative zonal component at the equatorial region which contributes to the flow through the first term in the RHS of eq. 1, but not through the second term where only nonzonal components of ξ appear. Then as ∇ H ξ points radially inwards, the corresponding flow circulation is eastward (in Fig. 4), the same direction as around the TC. The relative sign of the ξ streamfunction of mode 2 near the TC and at the equator appears in Figure 7 and the common direction of the zonal longitudinal flow appears in Figure 9. The two coupled spatial functions of mode 3 differ more than for modes 1 and 2 (Table 4). In particular, EOF 3 for flow guf m1 (Fig. 4 bottom left) is mainly localized under the Atlantic, a feature not seen in the EOF 3 that reconstructs flow COV -OBS (Fig. 4 bottom right). There are nonetheless some corresponding features in these EOF 3 charts: a cyclone under West of Iberia and another under East of the Caribbean Sea, an anticyclone under Eastern Canada and another under the Bering Sea. As to the expansion coefficient time series (PCs), notice that for a given mode the flow guf m1 and flow COV -OBS PCs for the SVD analysis are always closer to each other than the PC's of separated EOF/PC analysis of flow guf m1 and flow COV -OBS (see Figure 5). In addition, for mode 1 (respectively mode 2), the expansion coefficients obtained by SVD are closer to the one obtained by PCA of flow COV -OBS (respectively flow guf m1 ). This ability of SVD to bring the two flow descriptions closer together is also observed for spatial structures: the EOFs for the SVD analysis shown in Figure 4 are closer together than the corresponding EOFs for PCA analysis of individual flows (see Figure 2). This result is also apparent in Table 4, from comparison of values inside and outside square brackets for congruence (spatial similarity) and correlation (temporal resemblance) coefficients.

A quasi-periodicity is seen in the temporal variation of mode 2. Given the importance of this mode in terms of data variability explained and significance of its spatial features close to the TC and at the equator, we further characterize this temporal behaviour. To this end, all P C 2 s obtained from PCA applied to flow guf m1 and flow COV -OBS separately, both in a global grid and hemispherical subdomain grids, together with P C 2 s obtained from SVD of coupled flows, in a total of 9 time series, were fitted with a single sinusoid where amplitude, frequency and initial phase were free parameters. An average and a standard deviation of the estimated period T for the ensemble PC 2 s were produced. When using the whole 1840-1990 dataset we found T = 90.5 ± 4.6 yr and, when using the more recent dataset for for which errors in geomagnetic field models are smaller, we found T = 80.1 ± 2.7 yr.

Correlating PC i of flow guf m1 (resp. flow COV -OBS ) with the data matrix corresponding to flow guf m1 (resp. flow COV -OBS ) yields an homogeneous correlation map, whereas correlating it with the data matrix corresponding to flow COV -OBS (resp. flow guf m1 ) yields an heterogeneous correlation map (see Figure 6). Homogeneous correlation maps provide a representation of EOF i , while heterogeneous correlation maps are of special interest, since they reveal the correlated structures present in the two flows. Strong blue and red patches show regions of flow guf m1 or flow COV -OBS where the flow evolution is strongly correlated with a certain PC i (red for positive correlation, blue for anti-correlation). Note that the black solid lines, which separate regions of positive and negative correlation with PC i , are also the zero iso-ξ lines of EOF i .

The correlation maps of Figure 6 highlight a spiraling structure in mode 1, better seen in flow guf m1 than in flow COV -OBS and that was not obvious in Figure 4 because of disparate flow amplitudes. It is brought out by the normalization underlying the computation of correlation coefficients.

If we denote the mean flow column vector as x, with the same dimension as u i and v i vectors (the EOFs to retrieve flow guf m1 and flow COV -OBS , respectively), then we may write the mean flows as a combination of orthogonal SVD modes x guf m1,T = α T U T and x COV -OBS,T = β T V T . The column vectors α and β contain the expansion coefficients of the mean flows in the orthogonal basis of EOFs. Note that, because the mean flows x guf m1 and x COV -OBS have been removed prior to PCA or SVD analysis applied to matrices X, there is no orthogonality imposed between them and the variability modes. , respectively (see text). The last row shows values of u rms mean and standard deviation for the whole flows flow guf m1 and flow COV -OBS (first and second column). Note that they induce a SV signal that doesn't fit exactly the SV models, mainly because of the parameterization error (see appendix A). The percentage of SV root mean square explained by both global flows is on average 87%.

The cumulative importance of the three first SVD modes both in terms of u rms and the energy in the secular variation that they account for is gathered in Table 5. We can conclude that the compressed description of the flows using the mean and the three first variability modes, already accounts for about 90% of the total u rms . The percentage of SV explained by this simplified flow description is even greater, amounting to 95% or more. To see this, each mode added to the mean flow and previous variability modes, was made to interact with the corresponding geomagnetic field model for the whole 1840-1990 period and the relative SV energy was computed for each epoch. From Table 5 we can conclude that the variability in flow COV -OBS is higher than in flow guf m1 , since larger fractions of the whole flow u rms and induced SV are explained by variability modes. Finally, it is also shown in the last two columns the fraction of mean flow projected onto each one of the three modes. Mode 2, in particular, contributes significantly to both mean flows.

Estimates of significance of the correlation coefficients were computed as explained in appendix B. To build one single map that could concentrate information from the SVD analysis, we computed for mode i, with i = 1, 2, 3, the mean of the two PCs, PC guf m1 i and PC COV -OBS i since they are very close (see Figure 5). Then, the resulting time function P C i was correlated with flow guf m1 and with flow COV -OBS , giving two correlation maps r guf m1 i and r COV -OBS i , respectively. To condense information in one single map, a new parameter (r guf m1 i r COV -OBS i ) 1/2 (r guf m1 i /|r guf m1 i |) is introduced and its geographical distribution over the northern hemisphere is shown in Figure 7 using a color code. For each of the correlation coefficients r guf m1 i and r COV -OBS i , a corresponding p-value chart was computed (see sec. B). The p = 0.1 iso-contours for correlation of P C i (t) with flow guf m1 (dashed lines) and for correlation of P C i (t) with flow COV -OBS (solid lines) are plotted on Fig. 7. These contours enclose the regions with correlation coefficients r guf m1 i or r COV -OBS i statistically significant at level 90% or higher (i.e, the probability that these correlations are obtained by chance is no more than 10%), which means their variability is well explained by P C i (t) in both flow models. In deep blue (or red) regions, both flows evolve in time strongly correlated (or anticorrelated) with the corresponding P C i (t). Furthermore, if such patches are inside p = 0.1 iso-contours of both r guf m1 i and r COV -OBS i , we can conclude with a relatively high confidence (i.e., a 90% probability of success) that they represent spatial features of variability modes common to both flows, their temporal variation being represented by P C i (t). From Figure 7 top, we conclude that not only the main 3 vortices are significant features of the first mode, but also several smaller ones. Likewise, the zonal equatorial flow characterizing mode 2 is also a relevant feature in both flows, unlike the zonal flow component next to the tangent cylinder which is significant in flow guf m1 but not in flow COV -OBS . As explained above, the longitudinal flow component is in the same direction near the TC and in the equatorial region, in spite of different signs of ξ there (resulting in different signs of the correlation with P C 2 (t)). As for mode 3, the presence of correlated rolls under AH is significant, cyclonic around and 1970, anticyclonic around 1920 (see also Fig. 5). Their localization is however fuzzy as we can conclude from a relatively weak superposition of corresponding p-level curves from correlation with flow guf m1 and flow COV -OBS . Under the PH, this third mode shows more variability in flow COV -OBS than in flow guf m1 .

In order to visualize how and when these significant flow structures can be seen in the flow, Figure 8 shows the flow reconstructed from the mean and the two first SVD modes for flow guf m1 , at epochs 1860, 1900, 1930, 1970 and 1980. Referring only to the most significant features simultaneously spotted in both flows (see Fig. 7), the jet breakup can occur either because of an intensification of the cyclone centered at ∼ -150 • E longitude (mode 1), either because of an intensification of the high latitude cyclones cutting the TC under Asian continent, the most important centered at ∼ 60 • E longitude (mode 2). Negative values of PC and/or PC 2 tend to reinforce the jet, as can be seen in 1930.

Discussion

Mean flow

The mean flows obtained from our inversions of gufm1 and COV-OBS, shown on Fig.

are not a surprise, both displaying an eccentric gyre already described by [START_REF] Pais | Quasi-geostrophic flows responsible for the secular variation of the earth's magnetic field[END_REF]. There are however several differences between the two: flow COV -OBS needs a much weaker Pacific recirculation, and a lower overall rms value for the mean flow. As can be seen, the centennial timescale flow representing the mean shows strong latitudinal jets as main features, an anisotropy which is seen on different natural rotating flow systems, with or without magnetic fields. This fact was used by [START_REF] Schaeffer | On symmetry and anisotropy of Earth-core flows[END_REF] to propose an anisotropic regularization favoring zonation and leading to a mean flow not much different from that in Figure 1. The dominance of m = 1 azimuthal wavenumber in the mean flow suggests that an m = 1 forcing strongly influences the core flow. With direct numerical simulations of the core, [START_REF] Aubert | Bottom-up control of geomagnetic secular variation by the Earth's inner core[END_REF] advocate for a differential inner-core growth [START_REF] Monnereau | Lopsided Growth of Earth's Inner Core[END_REF][START_REF] Alboussière | Melting-induced stratification above the Earth's inner core due to convective translation[END_REF] to explain the jet eccentricity. In their model, the westward circulation is an effect of uppermost liquid core lagging behind the mantle which, due to gravitational coupling with the inner core is indirectly pushed eastward by thermochemical winds inside the tangent cylinder. However, the pattern observed here is strikingly close to the one found by [START_REF] Hori | Ancient dynamos of terrestrial planets more sensitive to core-mantle boundary heat flows[END_REF] (compare their Fig 4c with our Fig. 1) with a geodynamo simulation with internal sources and an heterogeneous (m = 1) heat flux imposed by the mantle. Thus the mean flows inverted from both COV-OBS and gufm1 could suggest a control of the convection by heterogeneous heat extraction from the mantle. They might indicate that more heat is extracted from the core in an hemisphere centered on the Western Atlantic, between -60 • and -90 • longitude, if we follow the findings of [START_REF] Hori | Ancient dynamos of terrestrial planets more sensitive to core-mantle boundary heat flows[END_REF].

First mode

The principal component analysis of both flow models show similar empirical orthogonal modes. The singular value decomposition of the coupled flows further pinpoints common features. The method is particularly convenient in making corresponding time functions (PCs) and spatial patterns (EOFs) that reconstruct flow guf m1 and flow COV -OBS more similar (Table 4). The empirical mode that carries the largest variability of the flow consists mainly of three big vortices (see Fig. 4): two anti-cyclones located at mid-latitudes around -60 • and +45 • longitude respectively; a cyclone located under the Eastern Pacific ocean around -150 • longitude, but which appears weaker on the flow inverted from COV-OBS. The significance analysis carried out on this mode further indicates that all the main three rolls (as well as a few weaker ones) are meaningful, the centers of the most significant regions being located at medium latitudes between 40 • and 50 • (see Fig. 7). Note also the cyclonic features close to the tangent cylinder under the North Atlantic ocean. The main effect of this empirical mode when added to the mean flow is to slightly reinforce the mean eccentric gyre before 1940, and later break it by making the three vortices largely dominant after 1960, and especially around 1970 (see Fig. 8).

Second mode

The second mode coming out of our analysis has the interesting property of carrying most of the angular momentum of our core flows, as shown in Fig. 9 (right) where the estimations for length of day (LOD) variations δT have been computed for our core flows and compared to LOD observations. The flow contributing to these variations is concentrated near the equator and the tangent cylinder (see Fig. 9, left). In addition to the zonal part, mode 2 also shows several vortices and especially two co-rotating ones located close to the tangent cylinder under Eastern and Western Russia. There are also vortices with opposite circulation, the most significant one being located under north of the Mediterranean basin, at mid latitudes (see Figure 7, region with 10 • E < long. < 20 • E and 30 • < lat. < 50 • N). During the period 1840-1990, this mode shows two main oscillations basically reinforcing the eccentric gyre around 1850 and 1930, and weakening it especially around 1900 with vortices and an opposite circulation around the tangent cylinder. The fact that most of the variation in zonal motion is correlated to these vortices might indicate either an excitation of a (zonal) eigenmode by the convecting flow, or the presence of slave secondary vortices excited by the interaction of a zonal eigenmode with the background magnetic field.

In this study, we show that these oscillations, have a mean period between 80 and 90 years (see Fig. 5). These rather long periods cannot be explained by torsional oscillations which are thought to have much shorter periods of about 6 years [START_REF] Gillet | Fast torsional waves and strong magnetic field within the Earth's core[END_REF], but could, in principle, be a match for (a combination of) MAC oscillations [START_REF] Braginsky | MAC-oscillations of the hidden ocean of the core[END_REF], as shown Figure 9: On the left, the zonal component of u φ for SVD mode 2, in km/yr, as a function of normalized distance to the rotation axis. On the right, the observed LOD variation (thick black line with squares), estimates from our SVD mode 2 assuming a columnar flow [thick blue and green lines, see Jault and Finlay, accepted, for the formula] or assuming the flow is restrained to a stratified layer at the top of the core of thickness 140 km (thin lines). recently by [START_REF] Buffett | Geomagnetic fluctuations reveal stable stratification at the top of the Earth's core[END_REF]. However, these MAC oscillations are restricted to a thin stratified layer at the top of the core. Fig. 9 (right) shows that the simple restriction of our QG core flow to a thin layer largely underestimates the LOD variations. Hence, with the hypothesis of MAC oscillations, the core flow cannot explain the LOD variations directly (at least not for the rms values of 10-20 km/yr constrained by the SV), while a QG model extending in the whole core nicely explain these [e.g. Gillet et al., 2010, Schaeffer and[START_REF] Schaeffer | On symmetry and anisotropy of Earth-core flows[END_REF]. It might be worth mentioning that the solar activity exhibits a period of 80 to 100 years, the Gleissberg cycle [e.g. Hathaway, 2010], which might provide an excitation for our mode 2. Periods of 96 yrs and 76 yrs associated to damped waves with parameters fitted to core surface flows have been reported by Zatman andBloxham [1997, 1998]. A core flow mode with a similar period has also been identified by [START_REF] Dickey | Leading modes of torsional oscillations within the Earth's core[END_REF] (85 yrs) and [START_REF] Buffett | Inversion of torsional oscillations for the structure and dynamics of Earth's core[END_REF] (86.3 yrs) using different methods. However, all those studies focused on the zonal longitudinal flow component inverted from gufm1 or a prior version, ufm1 [START_REF] Bloxham | Time-Dependent Mapping of the Magnetic Field at the Core-Mantle Boundary[END_REF]. In some cases, a shorter time span was analysed [e.g. Zatman andBloxham, 1997, 1998]. The standing oscillation reported in our study also involves nonzonal features, mainly close to the TC, and is detectable during the whole 1840 -1990 period. Furthermore, it combines information retrieved from two geomagnetic field models, gufm1 and COV-OBS. It must be emphasized that the time-dependence of mode 2 is not a nice sine wave, and that several frequencies are needed to reconstruct the signal. We obtain a dominant period of 80 to 90 years for the time-span 1840-1990 as well as for SVD applied to the shorter interval 1900-1990. However, if we restrict the SVD analysis from 1920 to 1990 the period decreases to ∼ 70 years and a 60 years sinusoid fitted to mode 2 gives an almost equally high R-squared value, in agreement with other studies [START_REF] Roberts | On the 60-year signal from the core[END_REF][START_REF] Buffett | Geomagnetic fluctuations reveal stable stratification at the top of the Earth's core[END_REF].

Flow reconstruction and other modes

The flow reconstructed from the first two modes and the mean flow is shown in Figure 8 and basically exhibits a very strong eccentric gyre around 1860 and 1930, which is destroyed by the appearance of several vortices around 1900 (especially two counter rotation vortices under Asia) and by three very strong and large vortices around 1970, basically splitting the unique global gyre into two circulation systems, a cyclonic one located under the Pacific ocean, and an anticyclonic one under Europe and the Atlantic ocean.

The third mode is less important in terms of the pseudo-streamfunction ξ variability explained, still it accounts for more than 10% either using PCA or SVD (Fig. 3). It's spatial pattern is also more complex showing a distribution of smaller-scale vortices, with often no clear correspondence between flow guf m1 and flow COV -OBS . The correlated vortices clustered under the Atlantic region show nonetheless some correspondence in both flows, significant only in reduced regions of Fig. 7.

We also checked higher order modes, which concentrate smaller fractions of the data variability. We looked in particular for a 6 year periodicity which has been reported by [START_REF] Gillet | Fast torsional waves and strong magnetic field within the Earth's core[END_REF]. Such periodicities can be glimpsed at modes 8 and 10 of PCA applied to COV-OBS for the time period 1960-2010. However, they never stand out above other relatively small periods that are also present, as is the case of a period of 9 years. Indeed, to isolate this component, [START_REF] Gillet | Fast torsional waves and strong magnetic field within the Earth's core[END_REF] need to extract the zonal flow of the ensemble average which is then band-pass filtered around 6 years, resulting in a signal several order of magnitudes smaller than the whole flow. Other tools should be employed to look at small signals of specific periods, as the EOF/PC analysis is not the adequate tool to isolate such a relatively low amplitude variability.

Conclusions

We use Principal Component Analysis tools applied to 'data' consisting of the pseudostreamfunction ξ evaluated on a regular spatial grid at the core surface outside the tangent cylinder. Our main goal is to identify some underlying structure in the flow and reconstruct the data from a linear combination of a small number of spatial patterns multiplied by timevarying coefficients. In this attempt for a simplification of observations, we first extracted the mean flow over the observatory era period (1840 to present): a large eccentric jet, flowing at lower latitudes under the Atlantic Hemisphere (AH) and at higher latitudes under the Pacific Hemisphere (PH) and a large cyclone under the PH, centered at medium-latitudes. The clear dichotomy revealed in this mean flow can be due to heterogeneous thermal forcing at the core-mantle boundary or at the solid core boundary [START_REF] Aubert | Detecting thermal boundary control in surface flows from numerical dynamos[END_REF], Aubert et al., 2013[START_REF] Hori | Ancient dynamos of terrestrial planets more sensitive to core-mantle boundary heat flows[END_REF].

We further identified three main circulation modes that account together for about 70% of the observed variability, 90% of the total flow u rms and 95% of the SV predicted by the total flow. Mode 1, consists of three large QG vortices indicative of large scale convection. It has an aperiodic variability during the inspected period, with a main boost around the 1970 epoch. At this time, the amplitude of mode 1 largely dominates the mean flow u rms and the three-roll flow breaks the large eccentric jet. The same dichotomy as in the mean flow reveals the dominance of anticyclonic circulation under the AH and cyclonic one under the PH.

A second mode was identified, which concentrates variations of core angular momentum, with a quasi-periodicity of 80 to 90 yrs. Periods close to 80 or 90 yrs have been found in magnetic field observatory data [START_REF] Jackson | Geomagnetic variation on decadal time scales: What can we learn from Empirical Mode Decomposition?[END_REF], geomagnetic field models [START_REF] Kang | Periodical Characteristics of the Geomagnetic Secular Variation and Decadal Length-of-Day Variation[END_REF] and inverted flows [START_REF] Zatman | Torsional oscillations and the magnetic field within the Earth's core[END_REF][START_REF] Kunsch | The jackknife and the bootstrap for general stationary observations[END_REF][START_REF] Dickey | Leading modes of torsional oscillations within the Earth's core[END_REF][START_REF] Buffett | Inversion of torsional oscillations for the structure and dynamics of Earth's core[END_REF]. Never before, however, have they been put forward without focusing the analysis on the zonal component of the flow and resorting to flows inverted from different geomagnetic field models. Here, this period emerges associated with an empirical orthogonal mode of the whole system, naturally associated to a large zonal flow but also, interestingly, to some significant small scale circulations mainly next to the tangent cylinder. These vortices evolve correlated with the equatorial flow and can also contribute to destroy or reinforce the large jet. The spatial structure of this mode has an important projection on the mean flow suggesting, if we believe on the presence of this flow component inside the core, that it can give some insight into the morphology of the sources responsible for accelerating the large eccentric jet. The physical existence of this mode needs further testing, however. The presence of a mode of 80-90 yrs period both in historical observatory data series and field models spanning centennial time periods seems confirmed. However, its absence from decade LOD historical data where a period of ∼ 60 yr is seen instead can sound suspicious, if we accept that decade LOD variations have their origin in angular momentum exchanges between the mantle and the core. A possible explanation would be that this flow mode results from a leakage of an external magnetic field signal (associated to the Gleissberg Sun-cycle) into internal geomagnetic field models and inverted core flows. Unfortunately, we are not in a position to deny (or confirm) this hypothesis based only on results from this study.

Mode 3 shows the presence of a few correlated rolls under AH but their localization is fuzzy due to the fact that they are retrieved differently in flow guf m1 and flow COV -OBS . This precludes the emergence of significant spatial features.

Significance tests were carried out with results that support the spatial stability of modes 1, 2 and 3 in general, and of specific spatial features present in these modes in particular. The PCA tools used in this study look for the directions of maximum variance in data space. The modes corresponding to projection along these directions will be of special interest in the case that higher variability modes correspond to interesting dynamics and lower ones correspond to noise. We can not guarantee from our results that this condition is fulfilled. However, we point out the similarities we found between certain identified modes (modes 1 and 2), significant both in terms of spatial features described and amount of variance explained, and other time/spatial features identified in previous studies of core dynamics.

This study gives important insight into the variability of the core flows on centennial time-scales. We have shown that the mean flow and three empirical modes can account for most of the secular variation of the geomagnetic field. Direct physical interpretation of the observed variability modes is not straightforward, although we have been able to link the first one to a dominant vortex pattern, and the second one to observed length-of-day variations. These modes may also guide future analytical and numerical studies on QG modes [e.g. [START_REF] Canet | Hydromagnetic quasi-geostrophic modes in rapidly rotating planetary cores[END_REF] and maybe help to constrain geometry of the internal magnetic field.

Finally, decomposition of quasi geostrophic flows inverted from geomagnetic field models into a few number of spatial structures and time dependent coefficients is of practical interest to serve as input for dynamical studies of the Earth's core interior where conditions closer to geomagnetic data are sought. An interesting question that this decomposition in empirical modes could help to answer is whether and how the core flow that we can reconstruct from magnetic models contributes to the generation of the Earth's magnetic field. We plan to address this question in a future study. observables and theoretical predictions and a term with quadratic forms on the flow model u is minimized [e. g Gubbins, 1983]:

Φ(m) = A(b) m -ḃ T C -1 e A(b) m -ḃ + λ m T C -1 m m . ( 13
)
m is the vector of poloidal (s m(c,s) ℓ

) and toroidal (t m(c,s) ℓ

) coefficients of the spherical harmonic (SH) expansion of the scalars S and T that determine u = R c ∇ H S -R c r ∧ ∇ H T ; ∇ H = ∇ -r∂/∂r is the horizontal nabla operator; b is the vector of SH time-dependent coefficients of the scalar potential V that determines the internal component of the geomagnetic field at the Earth's surface B = -∇V and ḃ is the vector of the corresponding first time derivatives; A is the interaction matrix or matrix of equations of condition, with elements that depend on the geomagnetic field model and on the Elsasser and Adams-Gaunt integrals [e.g. [START_REF] Whaler | Geomagnetic evidence for fluid upwelling at the core-mantle boundary[END_REF]:

ḃ = A(b) m + e . ( 14
)
The solution to this minimization problem is of the form

m = A T C -1 e A + λC -1 m A T C -1 e ḃ . ( 15
)
The formalism above is standard and can be found in different studies [see e.g. [START_REF] Pais | Length of day decade variations, torsional oscillations and inner core superrotation: evidence from recovered core surface zonal flows[END_REF]. We now concentrate on specific features of our inversions.

A.1 The error covariance matrix, C e

In the inversions for the present study, coefficients up to degree and order ℓ b = ℓ ḃ = 13 from model mod are used, for the deterministic part of the main field and its SV, where mod stands for gufm1 or COV-OBS. As to the dimension of vector m, it is that required to solve for flow coefficients up to degree and order ℓ m = 26, the maximum degree that can be constrained (even if only slightly) by the first 13 spherical harmonic degrees of the MF and of the SV.

The error e is treated as a Gaussian random variable with mean zero and covariance matrix C e , which can be condensed using the notation e ∼ N (0, C e ). We consider two contributions for C e : (i) C mod e is the a posteriori covariance matrix for the SV coefficients coming from the geomagnetic field model; (ii) C r e is the covariance matrix characterizing the spatial resolution or modeling errors. This second term is dependent on the flow solution m:

C e = C mod e + C r e (m) . (16)
For C mod e , we will be using the diagonal C guf m e matrix from [START_REF] Jackson | Time-dependency of tangentially geostrophic core surface motions[END_REF] and the dense C COV -OBS e matrix provided by [START_REF] Gillet | Stochastic modeling of the earth's magnetic field: Inversion for covariances over the observatory era[END_REF]. Previous studies have identified the modeling errors arising from an incomplete knowledge of the main field (C r e (m)) as having a higher contribution to the large scales of the SV than do observational errors [e.g. [START_REF] Eymin | On core surface flows inferred from satellite magnetic data[END_REF][START_REF] Pais | Quasi-geostrophic flows responsible for the secular variation of the earth's magnetic field[END_REF][START_REF] Gillet | Ensemble inversion of time-dependent core flow models[END_REF]. Inversions reported in recent papers do already take into account these errors, though different approaches have been followed, in particular the iterative approach in [START_REF] Pais | Quasi-geostrophic flows responsible for the secular variation of the earth's magnetic field[END_REF] and the stochastic approach in [START_REF] Gillet | Ensemble inversion of time-dependent core flow models[END_REF]. More recently another method has been proposed, called inverse geodynamo modeling, whereby the estimation of the effects of underparametrization is based on a statistical study of a numerical dynamo used as a prior model [START_REF] Aubert | Flow throughout the Earth's core inverted from geomagnetic observations and numerical dynamo models[END_REF][START_REF] Aubert | Earth's core internal dynamics 1840-2010 imaged by inverse geodynamo modelling[END_REF].

A law is required to prolongate the spectrum of the geomagnetic field to length scales that are unperceivable at the Earth's surface and above. We used the form R(ℓ, R c ) = W/(2ℓ + 1) for the Lowes-Mauersberger spectrum at the CMB, from [START_REF] Mcleod | Spatial and temporal power spectra of the geomagnetic field[END_REF]. This function was fitted to the geomagnetic spectrum computed from gufm-sat-E3 [START_REF] Finlay | Core surface magnetic field evolution 2000-2010[END_REF] for the 2005 epoch and using only 3 ≤ ℓ ≤ 13. The value W = 1.11 × 10 11 nT 2 /yr 2 was found. The small-scale magnetic field (14 ≤ ℓ ≤ 40) was then treated as a Gaussian random variable of zero mean and variance-covariance matrix C b S , i.e.

b S ∼ N (0, C b S ) , (17)
where C b S is a diagonal matrix with elements given by

σ 2 b S (ℓ) = W (ℓ + 1)(2ℓ + 1) 2 R c R E 2ℓ+4 , (18)
R E being the mean Earth radius.

To compute C e for a certain flow solution m, an ensemble of K modeling error vectors e r,k (m) = A(b S,k ) m are calculated, from an ensemble of K small-scale magnetic fields b S,k with statistical properties given by ( 17) and ( 18). From this sampling, estimates of the mean and covariance matrix for the modeling (or spatial resolution) error are computed according to standard formulas, the (non-diagonal) covariance matrix being denoted C r e (m). The updated error covariance matrix is then computed from eq. 16. In all our computations we used K = 80.

A.2 The flow regularization

The second term in the RHS of ( 13), also referred to as model regularization, incorporates the assumptions on the flow; its relative weight in the linear combination is controlled by the magnitude of the regularization parameter λ. Two important assumptions are considered in this study: (a) that most of the observed SV is due to a large scale flow interacting with the magnetic field and (b) that this flow is the surface expression of incompressible columnar convection inside the core. To impose condition (a) a ℓ 3 norm that minimizes flow gradients was used [START_REF] Gillet | Ensemble inversion of time-dependent core flow models[END_REF]. As to condition (b), it requires the flow at the CMB to be equatorially symmetric and verifying ∇ H • (u cos 2 θ) = 0 [e.g. [START_REF] Pais | Quasi-geostrophic flows responsible for the secular variation of the earth's magnetic field[END_REF][START_REF] Schaeffer | On symmetry and anisotropy of Earth-core flows[END_REF][START_REF] Amit | Differences between tangential geostrophy and columnar flow[END_REF]. This second condition is imposed using a very high λ multiplier, such that deviations are zero to the machine precision. The regularizing parameter multiplying the ℓ 3 norm, denoted λ R , is allowed to change from one epoch to the other, though, in order for the normalized misfit to be nearly 1 at each epoch. For each inversion the following condition is used to select an acceptable λ R [START_REF] Pais | Nonuniqueness of inverted core-mantle boundary flows and deviations from tangential geostrophy[END_REF]

: 0.8 < χ 2 N < 1.2 (19)
where

χ 2 = A(b) m -ḃ T C -1 e A(b) m -ḃ (20)
is the normalized misfit and N is the total number of SV coefficients.

A.3 The inversion algorithm

In this study, the iterative approach of [START_REF] Pais | Quasi-geostrophic flows responsible for the secular variation of the earth's magnetic field[END_REF] was used, whereby the effects of the small-scale magnetic field are computed for each step i of the flow calculation and the result is assimilated as a modeling error into the covariance matrix C i e . The procedure is repeated iteratively until a convergence is met, when either the modeling error or the flow solution no longer suffer significant changes. At this point, the found solution is consistent with the error covariance matrix C e used to invert it. Another iterative search is made with iteration index j, within each step i, in order to find the regularization parameter λ R that provides a flow solution m satisfying (19).

The whole procedure can be described in the following way, for each epoch. Each step i of the cycle involves three intermediate flows m i,0 , m i,1 and m i,2 in four intermediate steps within step i, denoted A) to D): A) m i,0 is used to compute a first estimate C i,0 e (m i,0 ) of the error covariance matrix, as explained in section A.1;

B) C i,0
e is used in eq. 15 to invert for m i,1 and from this a new error covariance matrix C i,1 e (m i,1 ) is computed as in previous item;

C) C i,1
e is used to invert for m i,2 and m i,2 is used to compute C i,2 e (m i,2 );

D) The flow solution resulting from step i is m i = (m i,1 + m i,2 )/2, which is used as m i+1,0 in the next step.

To monitor the convergence of flow solutions m i,1 and m i,2 and of covariance matrices C i,1 e (m i,1 ) and C i,2 e (m i,2 ) as i increases, they are compared using as diagnostic parameters the L 2 norm of the vector difference m i,2m i,1 and the Frobenius norm of the matrix difference C i,2 e (m i,2 ) -C i,1 e (m i,1 ). Each of the two flow inversions made at step i is done according to a second (internal) iterative cycle with index j, whereby the found solution must verify condition (19): j = 0, 1: Two trial solutions are first computed from (15), using λ R,0 and λ R,1 = λ R,0 + ∆λ 0 .

Computing χ 2 for each of these two inversions, makes possible to estimate the derivative of χ 2 with respect to λ. Then a Newton-Raphson type algorithm is applied to choose ∆λ 1 such that λ R,2 = λ R,1 + ∆λ 1 makes χ 2 approach the due condition (19).

1 < j < J: For each new value λ R,j = λ R,j-1 + ∆λ j-1 a new inversion using ( 15) is made and the corresponding χ 2 computed. A Newton-Raphson type algorithm is applied to choose ∆λ j . j = J: The cycle ends at iteration J, when λ R,J verifies condition (19). This final estimate gives λ R of the inversion.

In the first of all j-cycles, we used in our calculations λ R,0 = 1.0 × 10 6 and ∆λ 0 = λ R,0 /4, from previous tests. In all other cycles, the converged λ R from the previous cycle is used. Finally, the two particular i-cycle steps are: i = 0: The initial flow m 0,0 is computed from (15), using as error covariance matrix C e = C mod e + C r e , where only diagonal elements are considered for C r e , given by σ r e (ℓ) 2 = 36 exp(-ℓ) as in e.g. [START_REF] Schaeffer | On symmetry and anisotropy of Earth-core flows[END_REF].

i = I: Convergence is met for i = I if λ I,1 R = λ I,2
R AND χ I,1 (m I,1 ) = χ I,2 (m I,2 ) (see eq. 20), within a certain precision imposed by 19. These two conditions ensure that the inversion is well-converged, as confirmed by very small values of root squares of the L 2 norm of flow differences and of the Frobenius norm of error covariance matrices differences (10 -7 to 10 -8 and 10 -3 to 10 -5 , respectively.) Fig. 10 shows the Lowes-Mauersberger power spectrum at the Earth's surface of the SV of model COV-OBS, together with predictions from our inverted flows and different errors referred above: misfit errors, SV prior data errors from C COV -OBS e and modeling errors from C r e (m). For the two last cases, although the covariance matrices are dense, only diagonal elements are shown. Besides, rotationally invariant errors are calculated by averaging the diagonal elements (variances) within each degree ℓ. Along the time interval from 1840 to 1990 the prior data errors decrease as the quality and amount of data improve, but the modeling errors are always of the same order of magnitude with very similar spectra. It is worth noting that the prediction error for epoch 1850 is very similar to the corresponding error obtained by [START_REF] Aubert | Earth's core internal dynamics 1840-2010 imaged by inverse geodynamo modelling[END_REF], see his Figure 1], in spite of a more involved approach whereby the error covariance matrix to use is obtained from an ensemble of direct numerical simulations of the geodynamo.

When dealing with the two different geomagnetic field models gufm1 and COV-OBS, the main differences come from the fact that we used the diagonal C guf m e matrix from [START_REF] Jackson | Time-dependency of tangentially geostrophic core surface motions[END_REF] and the dense C COV -OBS e matrix provided by [START_REF] Gillet | Stochastic modeling of the earth's magnetic field: Inversion for covariances over the observatory era[END_REF]. The total matrices C e are always non-diagonal, contrary to most previous inversions, thus allowing to take into account the fact that errors in the SV coefficients are correlated.

For COV-OBS, this correlation comes in part from the geomagnetic field model computation, but the most important contribution comes from the parameterization errors iteratively computed during the inversion as explained in section A.1.

B Estimation of significances for identified modes

A test on the covariance matrix spectrum allows to assess the uncertainty of each eigenvalue of the covariance matrix C X (section 3.1) or each singular value of C X 1 X 2 (section 3.2), that reflects the fact that the sampling EOF's and PC's differ from the asymptotic ones (when N e → ∞). Following [START_REF] North | Sampling Errors in the Estimation of Empirical Orthogonal Functions[END_REF], the sampling error in the eigenvalue λ i of C X is, to first order,

δλ i ∼ λ i 2 N ( 21
)
where N is the number of independent realizations or the number of degrees of freedom, which we take equal to N e . Then, two consecutive eigenvalues λ i and λ j are considered degenerate if their difference ∆λ = λ iλ j is such that

∆λ ≤ λ i 2 N e . ( 22
)
This is known as North's rule of thumb. When two (or more) eigenvalues of C X are degenerate it means they cannot be clearly distinguished and the corresponding EOFs patterns cannot be considered independent [START_REF] Hannachi | Empirical orthogonal functions and related techniques in atmospheric science: A review[END_REF]. The same rule of thumb can be used to decide on the separation of consecutive singular values of the cross-covariance matrix when applying SVD [e.g. [START_REF] Venegas | Atmosphere-Ocean Coupled Variability in the South Atlantic[END_REF].

Temporal functions are compared using the correlation coefficient r [e.g. Von Storch and [START_REF] Storch | Analyses of Climate Variability-Applications of Statistical Techniques[END_REF], which measures the extent to which there is a linear relationship between the two functions. It always takes values in the interval [-1, 1], the two extreme values meaning there exists a perfect linear relationship between the two compared functions. For t 1 a column vector representing a certain time series and t 2 a column vector representing another time series, with the same dimension,

r(t 1 , t 2 ) = t 1 -t 1 T t 2 -t 2 t 1 -t 1 T t 1 -t 1 t 2 -t 2 T t 2 -t 2 (23)
where t 1 and t 2 are the average values of the two series over the time period considered.

To estimate the level of similarity between two different spatial structures, we use congruence coefficients, g [e.g. [START_REF] Harman | Modern Factor Analysis[END_REF]. If s 1 is a column vector representing a given spatial structure and s 2 a column vector for some other spatial structure, both with the same dimension, then

g(s 1 , s 2 ) = s T 1 s 2 s T 1 s 1 s T 2 s 2 (24)
Congruence coefficients have the same range as correlation coefficients, i.e. from -1 (negative similarity) to +1 (positive similarity). However, as noticed by [START_REF] Richman | Climatic Pattern Analysis of Three-and Seven-Day Summer Rainfall in the Central United States: Some Methodological Considerations and a Regionalization[END_REF] who uses this parameter in the context of PCA analysis, unlike the correlation coefficient the congruence coefficient does not subtract the mean value of the spatial patterns to compare. It therefore measures not only the pattern similarity but the magnitude similarity as well.

Finally, a Monte Carlo (MC) approach is used to estimate the significance of the correlation coefficients r for each grid point shown in correlation maps (see Figure 7). To this end, the two original time series to be correlated are shuffled, one at a time, using a technique called bootstrapping with moving blocks. The theoretical discussion of the method can be found in [START_REF] Kunsch | The jackknife and the bootstrap for general stationary observations[END_REF] and [START_REF] Liu | Moving blocks jackknife and bootstrap capture weak dependence[END_REF]. Here, we use it in a similar context as in several studies on coupled atmosphere-ocean variability [e.g. [START_REF] Wallace | Singular Value Decomposition of Wintertime Sea Surface Temperature and 500-mb Height Anomalies[END_REF][START_REF] Peng | The Coupled Patterns between Sea Level Pressure and Sea Surface Temperature in the Midlatitude North Atlantic[END_REF], Venegas et al., 1997]. The method is used for resampling and is applicable to weakly dependent stationary data that is, data which is nearly independent if far apart in time. There is nonetheless a correlation time to consider and that should be kept in resampled as in original data, i.e. the autocorrelation structure of the original series should be preserved. The reason is that shuffling the data in time randomly would increase the temporal degrees of freedom, giving in the end lower p-values for the statistical tests and improving artificially the significance of estimated parameters. So, in order to avoid this effect, the shuffling in time is applied not to every single data point randomly, but to a number of blocks containing the original data. The length of the blocks (b) depends on the autocorrelation structure of the original data. For a data set x i of total time-length N e , a number of blocks N eb + 1 are created so that the data from x 1 to x b will be in block 1, the data from x 2 to x b+1 will be in block 2, etc. After creating all the blocks, N e /b of these blocks are randomly selected to create a new pseudo-time series x ′ i , from which a certain quantity to test may be recalculated (if N e is not an integer multiple of b, x ′ i is truncated to have the same length as x i ). The length b of the blocks was chosen after inspection of the flow data autocorrelation function in each point of the spatial grid. Having verified that the time lag at which the autocorrelation function becomes zero or very closely zero does not exceed 55 yrs, this time interval was therefore used as the time-length of the moving blocks. For the significance tests on correlation maps produced from the SVD of coupled flow guf m1 and flow COV -OBS , each of the Monte-Carlo runs included two sub-runs. During the first sub-run the first of the series to correlate, (e.g. the P C i time series for SVD mode i, see section 3.2) was shuffled and the resulting artificial series was correlated in each grid point with the second series which is the observed time series from flow guf m1 or from flow COV -OBS . During the second sub-run the shuffling in time procedure was applied to the flow data matrix (flow guf m1 or flow COV -OBS , respectively), and the resulting resampled data series for each grid point were then correlated with P C i . The main goal is always to break the chronological order of one field relative to the other. For each MC run, the mean of the two artificial correlation coefficients for each grid point obtained from the two sub-runs was calculated. After a large number of MC runs (2 × 1000), the number of the correlation coefficients of magnitude higher than the original r divided by the total number of the runs is the p-value of this specific r. This p-value gives the probability to obtain this specific r just by chance. For instance, a p-value = 0.1 associated to a certain correlation coefficient, means that the probability to obtain such correlation coefficient for this specific pair of correlated time series just by chance is only 10%.

  mean QG-flow on the equatorial plane over the time period 1840 -1990, as seen from the North pole. For gufm1 (left) and COV-OBS (right). At the top, the pseudo streamfunction ξ/R c in units of km/yr; at the bottom, arrows visualize the flow.

Figure 2 :

 2 Figure 2: The three first EOFs resulting from PCA applied to flow guf m1 (left) and flow COV -OBS in the 1840-1990 period (right): mode 1 (top), mode 2 (middle) and mode 3 (bottom). Charted values are directly those from columns of matrices P (see section 3.1) and are normalized such that p T i p j = δ ij . Contours of ξ/R c streamfunction are used to visualize the flow, projected onto the equatorial plane, as seen from the North pole.

Figure 3 :

 3 Figure 3: Fraction of explained variablity f i for PCA applied to different datasets and for SVD of coupled flow guf m1 and flow COV -OBS (on bottom, right), with error bars with total length 2δf i , computed from eq. 21. Degenerate modes following North's criterion (eq. 22) are inside boxes .

  Fraction of variance, f i [%] (δf i [%]) (1840 -1990), PH 29.7 (3.7) 19.3 (2.4) 16.2 (2.0) 9.5 (1.2) 7.7 (1.0) COV-OBS (1840 -1990), AH 42.8 (5.5) 15.7 (2.0) 13.0 (1.7) 8.3 (1.1) 5.6 (0.7) gufm1 (1840 -1990), PH 40.3 (5.0) 23.4 (3.0) 11.5 (1.5) 7.4 (0.9) 5.3 (0.7) gufm1 (1840 -1990), AH 33.5 (4.3) 27.0 (3.4) 20.3 (2.6) 5.6 (0.7) 4.4 (0.6)

Figure 4 :

 4 Figure 4: The three first pairs of EOFs resulting from SVD of coupled fields applied to flow guf m1 and flow COV -OBS . For mode 1 (top), mode 2 (middle) and mode 3 (bottom), the EOFs allowing to reconstruct flow guf m1 (left) and flow COV -OBS (right). Values are directly those from columns of matrices U and V (see section 3.2) and are normalized such that u T i u j = v T i v j = δ ij . Contours of ξ/R c streamfunction are used to visualize the flow, projected onto the equatorial plane, as seen from the North pole.

Figure 5 :

 5 Figure 5: Temporal functions representing: the expansion coefficients of PCA modes that reconstruct flow guf m1 (thin solid line) and flow COV -OBS (thin line with empty circles); the expansion coefficients of SVD coupled modes that reconstruct flow guf m1 (thick solid line) and flow COV -OBS (thick line with filled circles).

Figure 6 :

 6 Figure 6: Homogeneous (two left columns, first for PC guf m1 i correlated with flow guf m1 , second for PC COV -OBS i correlated with flow COV -OBS ) and heterogeneous (two right columns, first for PC guf m1 i

Figure

  Figure 7: Maps for the parameters (r guf m1 i

Figure 8 :

 8 Figure 8: Contours of pseudo streamfunction ξ/R c in units of km/yr, for a reconstruction of flow guf m1 using the mean flow and the two first SVD modes. Shown epochs are, from top left to bottom right: 1860 (both P C 1 and P C 2 are weak, negative), 1900 (P C 2 strong, positive), 1930 (P C 2 strong, negative), 1970 (P C 1 very strong, positive) and 1980 (P C 1 and P C 2 strong, positive).In epochs of high positive values of PC 1 and/or PC 2 , the flow features in these two modes act to destroy the large eccentric jet, their amplitude being of the same order of magnitude as the mean flow. See also the movies[START_REF] Pais | COV-OBS (1840-2010): pseudostreamfunction; northern hemisphere[END_REF].

Table 2 :

 2 Subdomain stability of EOF/PCs for flow guf m1 and flow COV -OBS , when considering separately the Pacific and the Atlantic hemispheres. Percentage of variance explained by each mode and, in parenthesis, estimate of corresponding standard error, calculated as in Table1.

		gufm1	COV-OBS	gufm1	COV-OBS
	mode-1 0.99 / 0.95	0.87 / 0.77	0.94 / 0.95	0.98 / 0.98
	mode-2 0.88 / 0.77	0.62 / 0.65	0.93 / 0.94	0.93 / 0.96
	mode-3 0.85 / 0.81	0.96 / 0.90	0.96 / 0.97	0.92 / 0.90
	mode-4	0.25 / 0.30	0.68 / 0.79	0.72 / 0.70	0.39 / 0.34
	mode-5	0.25 / 0.32	0.87 / 0.85	0.73 / 0.84	0.21 / 0.28

|g(P H, global)| / |r(P H, global)| |g(AH, global)| / |r(AH, global)|

Table 3 :

 3 Subdomain stability of EOF/PCs for flow guf m1 and flow COV -OBS , when considering separately the Pacific and the Atlantic hemispheres. Absolute values of congruence coefficients between corresponding EOFs computed using the global grid and the grid under only the Pacific or the Atlantic hemispheres (|g(P H, global)| and |g(AH, global)|, respectively). Absolute values of correlation coefficients between corresponding PCs computed using the global grid and the grid under Pacific or Atlantic hemispheres (|r(P H, global)| and |r(AH, global)|, respectively). Bold is used for g and r values above 0.8, indicating high similarity.

	mode-1	f i [%] (δf i [%]) |g(EOF gufm1 i 33.9 (4.4)	, EOF COV -OBS i 0.80 [0.56] )| |r(P C gufm1 i	, P C COV -OBS i 0.96 [0.75] )|
	mode-2	19.0 (2.5)	0.74 [0.45]	0.93 [0.40]
	mode-3	16.4 (2.1)	0.62 [0.24]	0.94 [0.38]
	mode-4	8.4 (1.1)	0.61 [0.63]	0.92 [0.75]
	mode-5	7.1 (0.9)	0.47 [0.31]	0.91 [0.70]

Table 4 :

 4 SVD of coupled flows flow guf m1 and flow COV -OBS . In the first column, values for the percentage of cross-covariance explained by each coupled mode and, in parenthesis, estimates of corresponding standard error according to eq. 21. Second and third columns are for congruence (eq. 24) and correlation coefficients that compare corresponding EOFs and PCs, respectively, that reconstruct flow guf m1 and flow COV -OBS . In square brackets, values comparing corresponding EOFs and PCs from independent PCA applied to flow guf m1 and flow COV -OBS .

Table 5 :

 5 First and second columns show cumulative values of u rms,i = u

		u rms (km/yr)		Ḃi / Ḃ		mean flow projection
		gufm1 COV-OBS gufm1		COV-OBS		gufm1	COV-OBS
	mean flow mode-1 mode-2 mode-3 total flow	12.2 14.0 ± 3.2 12.2 ± 3.6 0.88 ± 0.19 0.86 ± 0.15 10.0 0.78 ± 0.16 0.75 ± 0.11 14.8 ± 3.1 13.1 ± 3.2 0.94 ± 0.20 0.91 ± 0.13 15.9 ± 4.0 14.6 ± 2.9 0.98 ± 0.17 0.95 ± 0.08 17.1 ± 5.1 16.7 ± 3.1 1.00 1.00	-0.17 0.35 0.22 -	-0.06 0.42 0.32 -
							i =
	(1/4πR 2 c ) u i • u i dS and standard deviation (std) values are displayed. Third and fourth columns show the cu-1/2 , starting from the mean flow over the 1840-1990 period; average
	mulative fraction of SV predicted, Ḃi / Ḃ where Ḃi = (1/4πR 2 c ) Ḃi • Ḃi dS comparing to the SV signal induced by the total flows, Ḃ ; average and std values are dis-1/2 , when
	played. Fifth and sixth columns show the normalized contribution of each orthogonal mode
	in the mean flow expansion, i.e., α 2 i /	Np k=1 α 2 k	1/2	and β 2 i /	Np k=1 β 2 k	1/2

Acknowledgments M. A. Pais is supported by FCT (PTDC/CTE-GIX/119967/2010) through the project COM-PETE (FCOMP-01-0124-FEDER-019978). A. L. Morozova is supported by a postdoc FCT scholarship (SFRH/BPD/74812/2010). M. A. Pais is grateful to the Grenoble University for funding her stay as an invited professor for two months. ISTerre is part of Labex OSUG@2020 (ANR10 LABX56). We thank two anonymous reviewers for helping us to improve this paper.

A Iterative inversion of gufm1 and COV-OBS

In this appendix we outline the procedure to invert the geomagnetic field models gufm1 in the time-period 1840-1990 and COV-OBS in 1840-2010, for a QG core flow with surface core expression u. The method is regularized weighted least squares inversion, whereby the objective function Φ(m) that expresses a linear combination of the discrepancy between