Determination of a high spatial resolution geopotential model using atomic clock comparisons

Recent technological advances in optical atomic clocks are opening new perspectives for the direct determination of geopotential differences between any two points at a centimeter-level accuracy in geoid height. However, so far detailed quantitative estimates of the possible improvement in geoid determination when adding such clock measurements to existing data are lacking. We present a first step in that direction with the aim and hope of triggering further work and efforts in this emerging field of chronometric geodesy and geophysics. We specifically focus on evaluating the contribution of this new kind of direct measurements in determining the geopotential at high spatial resolution (~ 10 km). We studied two test areas, both located in France and corresponding to a middle (Massif Central) and high (Alps) mountainous terrain. These regions are interesting because the gravitational field strength varies greatly from place to place at high spatial resolution due to the complex topography. Our method consists in first generating a synthetic high resolution geopotential map, then drawing synthetic measurement data (gravimetry and clock data) from it, and finally reconstructing the geopotential map from that data using least squares collocation. The quality of the reconstructed map is then assessed by comparing it to the original one used to generate the data. We show that adding only a few clock data points (less than 1 % of the gravimetry data) reduces the bias significantly and improves the standard deviation by a factor 3. The effect of the data coverage and data quality on the results is investigated, and the trade-off between the measurement noise level and the number of data points is discussed.


Introduction
Chronometry is the science of the measurement of time. As the time flow of clocks depends on the surrounding gravity field through the relativistic gravitational redshift predicted by Einstein (Landau and Lifshitz, 1975), chronometric geodesy considers the use of clocks to directly determine Earth's gravitational potential differences. Instead of using state-of-the-art Earth's gravitational field models to predict frequency shifts between distant clocks (Pavlis and Weiss (2003), ITOC project 1 ), the principle is to reverse the problem and ask ourselves whether the comparison of frequency shifts between distant clocks can improve our knowledge of Earth's gravity and geoid (Bjerhammar, 1985;Mai, 2013;Petit et al, 2014;Shen et al, 2016;Kopeikin et al, 2016). For example, two clocks with an accuracy of 10 −18 in terms of relative frequency shift would detect a 1-cm geoid height variation between them, corresponding to a geopotential variation ∆W of about 0.1 m 2 s −2 (for more details, see e.g. Delva and Lodewyck, 2013;Mai, 2013;Petit et al, 2014).
Until recently, the performances of optical clocks had not been sufficient to make applications in practice for the determination of Earth's gravity potential. However, ongoing quick developments of optical clocks are opening these possibilities. In 2010, Chou et al (2010) demonstrated the ability of the new generation of atomic clocks, based on optical transitions, to sense geoid height differences with a 30-cm level of accuracy. To date, the best of these instruments reach a stability of 1.6 × 10 −18 (NIST, RIKEN + Univ. Tokyo, Hinkley et al (2013)) after 7 hours of integration time. More recently, an accuracy of 2.1 × 10 −18 (JILA, Nicholson et al (2015)) has been obtained, equivalent to geopotential differences of 0.2 m 2 s −2 , or 2 cm on the geoid. Recently, Takano et al (2016) demonstrated the feasibility of cm-level chronometric geodesy. By connecting clocks separated by 15 km with a long telecom fiber, they found that the height difference between the distant clocks determined by the chronometric leveling (see Vermeer, 1983) was in agreement with the classical leveling measurement within the clocks' uncertainty of 5 cm. Other related work using optical fiber or coaxial cable time-frequency transfer can be found in (Shen, 2013;Shen and Shen, 2015).
Such results stress the question of what can we learn about Earth's gravity and mass sources using clocks, that we cannot easily derive from existing gravimetric data. Recent studies address this question; for example, Bondarescu et al (2012) discuss the value and future applicability of chronometric geodesy for direct geoid mapping on continents and joint gravity potential surveying to determine subsurface density anomalies. They find that a geoid perturbation caused by a 1.5 km radius sphere with 20 per cent density anomaly buried at 2 km depth in the Earth's crust is already detectable by atomic clocks with present-day accuracy. They also investigate other applications, for earthquake prediction and volcanic eruptions (Bondarescu et al, 2015b), or to monitor vertical surface motion changes due to magmatic, post-seismic, or tidal deformations (Bondarescu et al, 2015a,c).
Here we will consider the "static" or "long-term" component of Earth's gravity. Our knowledge of Earth's gravitational field is usually expressed through geopotential grids and models that integrate all available observations, globally or over an area of interest. These models are, however, not based on direct observations with the potential itself, which has to be reconstructed or extrapolated by integrating measurements of its derivatives. Yet, this quantity is needed in itself, like using a high-resolution geoid as a reference for height on land and dynamic topography over the oceans (Rummel and Teunissen, 1988;Rummel, 2002;Sansò and Venuti, 2002;Zhang et al, 2008;Rummel, 2012;Sansò and Sideris, 2013;. The potential is reconstructed with a centimetric accuracy at resolutions of the order of 100 km from GRACE and GOCE satellite data (Pail et al, 2011;Bruinsma et al, 2014), and integrated from near-surface gravimetry for the shorter spatial scales. As a result, the standard deviation (rms) of differences between geoid heights obtained from a global high-resolution model as EGM2008, and from a combination of GPS/leveling data, reaches up to 10 cm in areas well covered in surface data (Gruber, 2009). The uneven distribution of surface gravity data, especially in transitional zones (coasts, borders between different countries) and with important gaps in areas difficult to access, indeed limits the accuracy of the reconstruction when aiming at a centimetric level of precision. This is an important issue, as large gravity and geoid variations over a range of spatial scales are found in mountainous regions, and because a high accuracy on altitudes determination is crucial in coastal zones. Airborne gravity surveys are thus realized in such regions (Johnson, 2009;Douch et al, 2015); local clock-based geopotential determination could be another way to overcome these limitations.
In this context, here, we investigate to what extent clocks could contribute to fill the gap between the satellite and near-surface gravity spectral and spatial coverages in order to improve our knowledge of the geopotential and gravity field at all wavelengths. By nature, potential data are smoother and more sensitive to mass sources at large scales than gravity data, which are strongly influenced by local effects. Thus, they could naturally complement existing networks in sparsely covered places and even also contribute to point out possible systematic patterns of errors in the less recent gravity datasets. We address the question through test case examples of high-resolution geopotential reconstructions in areas with different characteristics, leading to different variabilities of the gravity field. We consider the Massif Central in France, marked by smooth, moderate altitude mountains and volcanic plateaus, and an Alps-Mediterranean zone, comprising high reliefs and a land/sea transition.
Throughout this work, we will treat clock measurements as direct determinations of the disturbing potential T (see below and Section 3 for details). We implicitly assume that the actual measurements are the potential differences between the clock location and some reference clock(s) within the area of interest. These measurements are obtained by comparing the two clocks over distances of up to a few 100 km. Currently two methods are available for such comparisons, fiber links (Lisdat et al, 2016) and free space optical links (Deschênes et al, 2016). The free space optical links are most promising for the applications considered here, but are presently still limited to short (few km) distances. However, projects for extending these methods based on airborne or satellite relays are on the way, but still require some effort in technology development.
The paper is organized as follows. In Section 2, we briefly summarize the method schematically. In Section 3, we describe the regions of interest and the construction of the high-resolution synthetic datasets used in our tests. In Section 4, we present the methodology to assess the contribution of new clock data in the potential recovery, in addition to ground gravity measurements. Numerical results are shown in Section 5. We finally discuss in Section 6 the influence of different parameters like the data noise level and coverage.

Method
The rapid progress of optical clocks performances opens new perspectives for their use in geodesy and geophysics. While they were until recently built only as stationary laboratory devices, several transportable optical clocks are currently under construction or test (see e.g. Bongs et al (2015); Origlia et al (2016); Vogt et al (2016). The technological step toward state-of-the art transportable optical clocks is likely to take place within the next decade. In parallel, in order to assess the capabilities of this upcoming technology, we chose an approach based on numerical simulation in order to investigate whether atomic clocks can improve the determination of the geopotential. Based on the consideration that ground optical clocks are more sensitive to the longer wavelengths of the gravitational field around them than gravity data, our method is adapted to the determination of the geopotential at regional scales. In figure 1 a scheme of the method used in this paper is shown: 1. In the first step, we generate a high spatial resolution grid of the gravity disturbance δ g and the disturbing potential T , considered as our reference solutions. This is done using a state of the art geopotential model (EIGEN-6C4), and by removing low and high frequencies. It is described in details in Section 3; 2. In the second step, we generate synthetic measurements δ g and T from a realistic spatial distribution, then we add generated random noise representative of the measurement noise. This is described in details in Section 4; 3. In a third step, we estimate the disturbing potential T from the synthetic measurements δ g and/or T on a regular grid thanks to Least-Square Collocation (LSC) method. Interpolating spatial data is realized by making an assumption on the a priori gravity field regularity on the target area, as described in Section 5. This prior is expressed by the covariance function of the gravity potential and its derivatives. It allows to predict the disturbing potential on the output grid from the observations using the signal correlations between the data points, and with the estimated potential.
4. Finally, we evaluate the potential recovery quality for different data distribution sets, noise levels, and types of data, by comparing the statistics of the residuals δ between the estimated values T and the reference model T .
Let us underline that in this work, we use synthetic potential data while a network of clocks would give access to potential differences between the clocks. We indeed assume that the clocks-based potential differences have been connected to one or a few reference points, without introducing additional biases larger than the assumed clock uncertainties. Note that these reference points are absolute potential points determined by other methods (GNSS/geoid for example).
Step 1: Build synthetic field model Step 2: Select data distribution and add noise Step 3: Make an assumption on the a priori gravity field and estimate a potential model In this differential method, significant residuals δ (higher than the machine precision) can have several origins, depending on the parameters of the simulation that can be varied: 1. The modeled instrumental noise added to the reference model at step 2. This noise can be changed in order to determine, for instance, whether it is better to reduce gravimetry noise by one order of magnitude, rather than using clock measurements; 2. The data distribution chosen in step 2. This is useful to check for instance the effect of the number of clock measurements on the residuals or to find an optimal coverage for the clock measurements; 3. The potential estimation error, due to the intrinsic imperfection of the covariance model chosen for the geopotential. In our case, this is due to the low-frequency content of the covariance function chosen for the Least-Square Collocation method (see Section 5).
All these sources of errors are somewhat entangled with one another, such that a careful analysis must be done when varying the parameters of the simulation. This is discussed in details in Section 6.
3 Regions of interest and synthetic gravity field reference models

Gravity data and distribution
Our study focuses on two different areas in France. The first region is the Massif Central located between 43 • to 47 • N and 1 • to 5 • E, and consists of plateaus and low mountain range, see Figure 2. The second target area, much more hilly and mountainous, is the French Alps with a portion of the Mediterranean Sea located at the limit of different countries and bounded by 42 • to 47 • N and 4.5 • to 9 • E, see Figure 3. Topography is obtained from the 30 m digital elevation model over France by IGN, completed with Smith and Sandwell (1997) bathymetry and SRTM data.  Available surface gravity data in these areas, from the BGI (International Gravimetric Bureau), are shown in Figures 2b-3b. Note that the BGI gravity data values are not used in this study, but only their spatial distribution in order to generate realistic distribution in the synthetic tests. In these figures, it is shown that the gravity data are sparsely distributed: the plain is densely surveyed while the mountainous regions are poorly covered because they are mostly inaccessible by the conventional gravity survey. The range of free-air gravity anomalies (see e.g Moritz, 1980;Sansò and Sideris, 2013) which are quite large reflects the complex structure of the gravity field in these regions, which means that the gravitational field strength varies greatly from place to place at high-resolution. The scarcity of gravity data in the hilly regions is thus a major limitation in deriving accurate high-resolution geopotential model.

High resolution synthetic data
Here, we present the way to simulate our synthetic gravity disturbances δ g and disturbing potentials T by subtracting the gravity field long and short wavelengths influence of a high-resolution global geopotential model.
The generation of the synthetic data δ g and T at the Earth's topographic surface was carried out, in ellipsoidal approximation, with the Fortran program GEOPOT 2 (Smith, 1998) of the National Geodetic Survey (NGS). This program allows to compute gravity field related quantities at given locations using a geopotential model and  additional informations such as parameters of the ellipsoidal normal field, tide system. The ellipsoidal normal field is defined by the parameters of the geodetic reference system GRS80 (Moritz, 1984). As input, we used the static global gravity field model EIGEN-6C4 (Förste et al, 2014). It is a combined model up to degree and order (d/o) 2190 containing satellite, altimetry, terrestrial gravity and elevation data. By using the spherical harmonics (SH) coefficients up to d/o 2000, it allows us to map gravity variations down to 10 km resolution. Thus, these synthetic data do not represent the full geoid signal. The choice is motivated by the fact that at a centimetric level of accuracy, we expect large benefit from clocks at wavelengths ≥ 10 km.
Our objective is to study how clocks can advance knowledge of the geoid beyond the resolution of the satellites. In a first step, as illustrated in Figure 4, the long wavelengths of the gravity field covered by the satellites and longer than the extent of the local area are completely removed up to the degree n cut = 100 (200 km resolution). This data reduction is necessary for the determination of the local covariance function in order to have centered data, or close to zero, as detailed in Knudsen (1987Knudsen ( , 1988. Between degree 101 and 583, the gravity field is progressively filtered using 3 Poisson wavelets spectra (Holschneider et al, 2003), while its full content is preserved above degree 583. In this way we realize a smooth transition between the wavelengths covered by the satellites and those constrained from the surface data. To subtract the terrain effects included in EIGEN-6C4, we used the topographic potential model dV_ELL_RET2012 (Claessens and Hirt, 2013) truncated at d/o 2000. Complete up to d/o 2160, this model provides in ellipsoidal approximation the gravitational attraction due to the topographic masses anywhere on the Earth's surface. The results of this data reduction yields to the reference fields δ g and T for both regions, shown in Figures 5-6.
The Figures 5-6 show the different characteristics of the residual field in these two regions. The residual anomalies have smaller amplitudes in the Massif Central area when compared to the Alps. In addition, the presence of high mountains on part of the latter zone results in an important spatial heterogeneity of the residual gravity anomalies, with large signals also at intermediate resolutions.

Dataset selection and synthetic noise
Gravimetric location points selection. Our goal is to reproduce a realistic spatial distribution of the gravity points. The BGI gravity data sets contain hundreds of thousands points for the target regions (see Figure 2b-3b).
In order to reduce the size of the problem and make it numerically more tractable, we build a distribution with no more than several thousand points from the original one. Starting from the spatial distribution of the BGI gravity data sets, a grid δ g of N cells is built with a regular step of about 6.5 km. Each cell contains n i points with i = {1, 2, . . . , N}. These n i points are replaced by one point which location is given by the geometric barycenter of the n i points, in the case that n i > 0. If n i = 0 then there is no point in the cell i. Figures 7 show the new distributions of gravimetric data for the Massif Central and the Alps regions; they have, respectively, 4374 and 4959 location points. These new spatial distributions reflect the initial BGI gravity data distribution but are be more homogeneous. They will be used in what follows.
(a) Massif Central: 4374 gravity data and 33 potential data.
(b) Alps: 4959 gravity data and 32 potential data. Chronometric location points selection. We choose to put clock measurements only where existing land gravity data are located. Indeed, these data mainly follow the roads and valleys which could be accessible for a clock comparison. Then, we use a simple geometric approach in order to put clock measurements in regions where the gravity data coverage is poor. Since the potential varies smoothly compared to the gravity field, a clock measurement is affected by masses at a larger distance than in the case of a gravimetric measurement. For that reason, a clock point will be able to constrain longer wavelengths of the geopotential than a gravimetric point. This is particularly interesting in areas poorly surveyed by gravity measurement networks. Finally, in order to avoid having clocks too close to each other, we define a minimal distance d between them. We chose d greater than the correlation length of the gravity covariance function (in this work λ ∼ 20 km, see Table 1).
Here we give more details about our algorithm to select the clock locations: 1. First, we initialize the clock locations on the nodes of a regular grid T with a fixed interval d. This grid is included in the target region at a setback distance of about 30 km from each edge (outside possible boundary effects); 2. Secondly, we change the positions of each clock point to the position of the nearest gravity point from the grid δ g, located in cell i (see the previous paragraph); in cell i are located n i points of the initial BGI gravity data distribution; 3. Finally, we remove all the clock points located in cells where n i > n max . This is a simple way to keep only the clock points located in areas with few gravimetric measurements.
This method allows to simulate different realistic clock measurement coverages by changing the values of d and n max . The number of clock measurements increases when the distance d decreases or when the threshold n max increases, and vice versa. It is also possible to obtain different spatial distributions but the same number of clock measurements for different sets of d and n max . In Figure 7, we propose an example of clock coverage used hereafter for both target regions with 32 and 33 clock locations, respectively, in the Massif Central and the Alps, corresponding to ∼ 0.7 percent of the gravity data coverage. For the chosen distributions, the value of d is about 60 km and n max = 15.
Synthetic measurements simulation. For each data point, the synthetic values of δ g and T are computed by applying the data reduction presented in Section 3.2. It is important to note that the location points of the simulated data T are not necessarily at the same place than the estimated data T .
A Gaussian white noise model is used to simulate the instrumental noise of the measurements. We chose, for the main tests in the next section, a standard deviation σ δ g = 1 mGal for the gravity data and σ T = 0.1 m 2 /s 2 for the potential data. In terms of geoid height, the latter noise level is equivalent to 1 cm. Other tests with different noise levels are discussed in Section 6.

Numerical results
In this section, we present our numerical results showing the contribution of clock data in regional recovery of the geopotential from realistic data points distribution in the Massif Central and the Alps. The reconstruction of the disturbing potential is realized from the synthetic measurements δ g and T , and by applying the Least-Squares Collocation (LSC) method.
Planar Least-Squares Collocation. The LSC method, described in Moritz (1972Moritz ( , 1980, is a suitable tool in geodesy to combine heterogeneous data sets in gravity field modelling. Assuming that the measured values are linear functionals of the disturbing potential T , this approach allows us to estimate any gravity field parameter based on T from many types of observables. Consider l = [ l T , l δ g ] = l k a data vector composed by p data T and q data δ g, affected by measurement errors ε k , with k = {1, 2, . . . , p + q}. The estimation of the disturbing potential T P at point P from the data l can be performed with the relation T P = C T P ,l · C −1 n,n · l (5.1) C n,n = C l,l + ω C ε,ε (5.2) with C l,l the covariance matrix of the measurement vector l, C ε,ε the covariance matrix of the noise, C T P ,l the cross covariance matrix between the estimated signal T P and the data l, and ω the Tikhonov regularization factor (Neyman, 1979), also called weight factor. In practice, the data l are synthesized as described in Sections 3 and 4. Therefore, the measurement noise is known to be a Gaussian white noise. Noise and signal (errorless part of l k ) are assumed to be uncorrelated, and the covariance matrix of the noise can be written as with I n the identity matrix of size n. Because C l,l can be very ill-conditioned, the matrix (5.3) plays an important role in its regularization before inversion, since positive constant values are added to the elements of its main diagonal. To avoid any iterative process to find an optimum value of ω in case where this matrix C l,l is not definite positive, we chose to fix the weight factor ω = 1 and to apply a singular value decomposition (SVD) to pseudo-inverse the matrix. As shown in (Rummel et al, 1979), these two approaches are similar.
Estimation of the covariance function. Implementation of the collocation method requires to compute the covariance matrices C T P ,l and C l,l . This step has been carried out using a logarithmic spatial covariance function from (Forsberg, 1987), see Appendix A. This stationary and isotropic model is well-adapted to our analysis.
Indeed, it provides the auto-covariances (ACF) and cross-covariances (CCF) of the disturbing potential T and its derivatives in 3 dimensions with simple closed-form expressions.
The spatial correlations of the gravity field are analyzed with the program GPFIT (Forsberg and Tscherning, 2008). The variance C 0 is directly computed from the gravity data on the target area, and the parameters α and β (see Appendix A) are estimated by fitting the a priori covariance function to the empirical ACF of the gravity disturbances δ g.
Results of the optimal regression analysis for both regions are given in Figure 8 and Table 1. The estimated covariance models reflect the different characteristics of the gravity signals in the two areas and the data sampling, which is less dense in high relief areas. Finally, the gravity anomaly covariances show similar correlation lengths, with a larger variance for the case of the Alps; their shapes, however, slightly differ, with a broader spectral coverage for the Alps.    Table 1: Estimation of the auto covariance function parameters on the gravity data δ g using the logarithmic model from Forsberg (1987) with, µ the mean, C 0 the variance, α and β respectively a shallow and a compensating depth parameter. Here, λ is the correlation length defined as the distance at which the covariance is half of the variance.
Knowing the parameter values of the covariance model, we can now estimate the potential anywhere on the Earth's surface.
Contribution of clocks. The contribution of clock data in the potential recovery is evaluated by comparing the residuals of two solutions to the reference potential on a regular grid interval of 10 km. The first solution corresponds to the errors between the estimated potential model computed solely from gravity data and the potential reference model, while the second solution uses combined gravimetric and clock data. To avoid boundary effects in the estimated potential recovery, a grid edge cutoff of 30 km has been removed in the solutions.
For the Massif Central region, the disturbing potential is estimated with a bias µ T ≈ 0.041 m 2 s −2 (4.1 mm) and a rms σ T ≈ 0.25 m 2 s −2 (2.5 cm) using only the 4374 gravimetric data, see Figure 9a. When we now reconstruct T by adding the 33 potential measurements to the gravimetric measurements, the bias is improved by one order of magnitude (µ T ≈ −0.002 m 2 s −2 or −0.2 mm) and the standard deviation by a factor 3 (σ T ≈ 0.07 m 2 s −2 or 7 mm), see Figure 9b.
For the Alps, Figure 10, the potential is estimated with a bias µ T ≈ 0.23 m 2 s −2 (2.3 cm) and a standard deviation σ T ≈ 0.39 m 2 s −2 (3.9 cm) using only the 4959 gravimetric data. When adding the 32 potential measurements, we note that the bias is improved by a factor 4 (µ T ≈ −0.069 m 2 s −2 or −6.9 mm) and the standard deviation by a factor 2 (σ T ≈ 0.18 m 2 s −2 or 1.8 cm).
It can be noticed that the residuals in both areas differ. This results from the covariance function that is less well modeled when the data survey has large spatial gaps. It should also be stressed that a trend appears in the reconstructed potential with respect to the original one when no clock data are added in both regions. This effect is discussed in Section 6.

Discussion
Effect of the number of clock measurements. Figure 11 shows the influence of the number of clock data in the potential recovery, and thereby of their spatial distribution density. We vary the number and distribution of clock data by changing the mesh grid size d, which represents the minimum distance between clock data points (see Section 4). The particular cases shown in detail in Section 5 are included. We characterize the performance of the potential reconstruction by the standard deviation and mean of the differences between the original potential on the regular grid and the reconstructed one. When increasing the density of the clock network, the standard deviation of the differences tends toward the centimeter level, for the Massif Central case, and the bias can be reduced by up to 2 orders of magnitude. Note that we have not optimized the clock locations such as to maximize the improvement in potential recovery. The chosen locations are simply based on a minimum distance and a maximum coverage of gravity data (c.f. Section 4). An optimization of clock locations would likely lead to further improvement, but is beyond the scope of this work and will be the subject of future studies.
Moreover, the results indicate that it is not necessary to have a large number of clock data to improve the reconstruction of the potential. We can see that only a few tens of clock data, i.e. less than 1 percent of the gravity data coverage, are sufficient to obtain centimeter level standard deviations and large improvements in the bias. When continuing to increase the number of clock data the standard deviation curve seems to flatten at the cm level.
(b) Alps area. Figure 11: Performance of the potential reconstruction (expressed by the standard deviations and mean of differences between the original potential on the regular grid and the reconstructed one) wrt the number of clocks. In green: number of clock data in terms of percentage of δ g data.
Effect of the number of gravity measurements. We have performed numerical tests in order to study the influence of the density of gravity measurements on the reconstructed disturbing potential, with or without clocks. We take the case of the Massif Central region, and set-up simulations where the clock coverage is fixed (either no clocks, or 38 clocks at fixed locations where we also have gravity data). Then, we progressively increase the spatial resolution of the gravity data, from 91 to 6889 points, and evaluate as before the quality of the potential reconstruction with or without clocks. Here, in contrast with the tests presented in the previous section, the gravity points are randomly generated from a complete 5-km step grid. Figure 12 shows the results of these tests. If we compare the rms values between configurations where we add clocks or not, we observe that the behavior of the results is globally similar, and improved with clocks. The interpolation error due to a too low resolution of the gravity data with respect to scales of the field variations predominates when we have less than ∼1500 gravity measurements, leading to large rms values even with clocks. Above this number, the large-scale reconstruction errors significantly contribute to the rms of residuals, explaining that the rms further decreases only when clocks are added. Looking at the bias between the reconstructed and original potential, we can see that it is poorly dependent on the number of gravity data in the tests without clocks. It probably reflects the fact that these data are more sensitive to the smaller scale components of the gravity potential. When we add clocks, the improvement on the bias is always important, which is consistent with the fact that the higher sensitivity of clocks to the longer wavelengths of the field reduces significantly the trend from the modelling error.  Covariance function consistency. In Figures 9a and 10a, a trend appears in the residuals, but disappears when gravimetric and clock data are combined. This is due to the fact that the covariance function does not have the same spectral coverage as the data generated from the gravity field model EIGEN-6C4. Indeed, the covariance function contains low frequencies while we have removed them for the synthetic data. Therefore, some low frequency content is present in the recovered potential. Whilst the issue could be avoided by using a covariance parametric model from which we can remove the low frequency content in a perfectly consistent way with the data generation (e.g. a closed-form Tscherning-Rapp model (Tscherning and Rapp, 1974;Tscherning, 1976)), it is not obvious that the corresponding results would allows realistic conclusions. Indeed, the spectral content of real surface observations, after removal of lower frequencies from a global spherical harmonics model, may still retain some unknown low frequencies. As consequence, it is not obvious to match to that of a single covariance function, while perfect consistency can only be achieved from synthetic data. We chose to keep this mismatch, thereby investigating the interest of clocks for high-resolution geopotential determination when our prior knowledge on the surface data signal and noise components is not perfect. More detailed studies on this issue are considered beyond the scope of our paper, which presents a first step to quantify the possible use of clock measurements in potential recovery.
Influence of the measurement noise. We have also investigated the effect of the noise levels applied to the synthetic data, see Tables 2-3, by using various standard deviations to simulate white noise of the measurements: σ T = {1, 0.1} m 2 s −2 for the clock measurements and σ δ g = {1, 0.1, 0.01} mGal for the gravimetric measurements. These results were obtained for the same conditions as in Section 5, i.e. 33 (resp. 32) clock data points and 4374 (resp. 4959) gravity data points for the Massif Central (resp. Alps).
We can see that adding clocks improves the potential recovery (smaller standard deviation σ and bias µ of the residuals) for both regions and whatever the noise of the gravimetric or clock measurements.
We observe that decreasing the noise of the gravity data by up to two orders of magnitude only improves the standard deviation of the residuals σ of the recovered potential by comparatively small amounts (less than a factor 2). This is probably due to the fact that the covariance function does not reflect the gravity field correctly in these regions, combined with a limited data coverage. Note that the low frequency content in the covariance function (see above) is unlikely to be the main cause here, as the comparatively small reduction of σ is also observed when clocks are present in spite of the fact that they remove the low frequency trend (c.f. figures 9b and 10b).
When adding clocks the standard deviations are decreased by up to a factor 3.7 with low clock noise (0.1 m 2 s −2 or 1 cm) and a factor 1.5 with higher clock noise (1 m 2 s −2 or 10 cm). The effect is stronger in the Massif Central region than in the Alps. We attribute this again to the mismatch between the covariance function and the complex structure of the gravity field, which is larger in the Alps.
Aliasing of the very high-resolution components. We have studied the aliasing of gravity variations at scales shorter than 10 km spatial resolution, that would be present in real data but under-sampled by the finite spatial density of the surveys. Errors in the topographic corrections may reach a few mGal for DTM (Digital terrain model) sampled at hundreds of meters resolution (Tziavos et al, 2009), while local geological sources may lead to gravity signals up to ∼ 10 mGal (Yale et al, 1998;Bondarescu et al, 2012;Castaldo et al, 2014). Furthermore, we have analyzed the Bouguer gravity anomalies from the BGI database along profiles in the Massif Central and the Alps, and found, after smoothing the profiles at 10 km resolution, high-resolution components with rms ∼ 1 mGal in the Massif Central, and ∼ 3 mGal in the Alps. An order of magnitude of the corresponding geoid variations can be derived by assuming that the gravity signals at a given spatial scale are created by a point mass at the corresponding depth. We find that a 5 km width, 5 mGal (resp. 10 mGal) gravity anomaly corresponds to a 1.3 cm (resp. 2.6 cm) geoid variation, above the centimetric level indeed.
We simulate these previously neglected signals beyond 10 km resolution by increasing the noise level on the gravity data in our tests, up to 5 mGal in the Massif Central, and 10 mGal in the Alps. Note that these rms values are large with respect to the observed high-resolution variabilities in the data. As previously, numerical simulations are performed with and without adding clocks, and the results are presented in the first column of Table 2-3. We can see that decreasing the accuracy of the gravimetric measurements increases the residuals as compared to the previous solutions. This is due to the fact that the signal-to-noise ratio decreases, degrading the covariance function modelling. However, our previous conclusions on the benefit of clocks remain the same, even in the presence of significant signals at the shortest spatial scales.

Conclusions
Optical clocks provide a tool to measure directly the potential differences and determine the geopotential at high spatial resolution. We have shown that the recovery of the potential from gravity and clock data with the LSC method can improve the determination of geopotential at high spatial resolution, beyond what is available from satellites. Compared to a solution that does not use the clock data, the standard deviation of the disturbing potential reconstruction can be improved by a factor 3, and the bias can be reduced by up to 2 orders of magnitude with only a few tens of clock data. This demonstrates the benefit of this new potential geodetic observable, which could be put in practice in the medium term when the first transportable optical clocks and appropriate time transfer methods will be developed (see Bongs et al, 2015;Lisdat et al, 2016;Deschênes et al, 2016;Vogt et al, 2016). Since clocks are sensitive to low frequencies of the gravity field, this method is particularly well-adapted in hilly and mountainous regions for which the gravity coverage is more sparsely distributed, allowing to fill areas not covered by the classical geodetic observables (gravimetric measurements). Additionally, adding new observables helps to reduce the modelling errors, e.g. coming from a mismatch between the covariance function used and the real gravity field.
In the same way, GPS and leveling data have been used, in combination with gravity data, to derive highresolution gravimetric geoids (Kotsakis and Sideris, 1999;Duquenne, 1999;Denker et al, 2000;Duquenne et al, 2005;Nahavandchi and Soltanpour, 2006). Using clocks is, however, different from performing GPS and leveling measurements. They provide an information of similar nature as the gravity data, in contrast with these geometric observations. The latter are affected by different sources of errors (e.g. Duquenne, 1998;Marti et al, 2001), and quite expensive in the case of leveling campaigns. We can expect that clocks could help identify and reduce errors in the gravity and GPS/leveling through their joint analysis for geopotential determination. Beyond the application considered in this work, the clocks can also contribute to the unification of height systems realizations (Shen et al, 2011;Denker, 2013;Shen et al, 2016;Kopeikin et al, 2016;Takano et al, 2016), connecting distant points to a high-resolution reference potential network.
To our knowledge, this is the first detailed quantitative study of the improvement in field determination that can be expected from chronometric geodesy observables. It provides first estimates and paves the way for future more detailed and in depth works in this promising new field.
To overcome some limitations in the a priori model, as discussed in the previous section, we intend in a forthcoming work to investigate in more details the imperfections of the covariance function model. Moreover, as the gravity field is in reality non-stationary in mountainous areas or near the coast, some numerical tests with non-stationary covariance functions will be conducted. Another promising source of improvement could be the optimization of the positioning of the clock data. For example, the correlation lengths and the variations of the gravity field could be used as constraints. A genetic algorithm could also be considered to solve this location problem. Finally, it will be interesting to focus on the improvement of the potential recovery quality by combining other types of observables such as leveling data and gradiometric measurements. As knowledge of the geopotential provides access to height differences, this could be a way to estimate errors of the GNSS technique for the vertical positioning, or contribute to regional height systems unification.