Geodetic VLBI with an artificial radio source on the Moon: a simulation study

We perform extensive simulations in order to assess the accuracy with which the position of a radio transmitter on the surface of the Moon can be determined by geodetic VLBI. We study how the quality and quantity of geodetic VLBI observations influence these position estimates and investigate how observations of such near-field objects affect classical geodetic parameters like VLBI station coordinates and Earth rotation parameters. Our studies are based on today’s global geodetic VLBI schedules as well as on those designed for the next-generation geodetic VLBI system. We use Monte Carlo simulations including realistic stochastic models of troposphere, station clocks, and observational noise. Our results indicate that it is possible to position a radio transmitter on the Moon using today’s geodetic VLBI with a two-dimensional horizontal accuracy of better than one meter. Moreover, we show that the next-generation geodetic VLBI has the potential to improve the two-dimensional accuracy to better than 5 cm. Thus, our results lay the base for novel observing concepts to improve both lunar research and geodetic VLBI.


Introduction
Very Long Baseline Interferometry (VLBI) has a long tradition of serving the astronomical, astrometric and geodetic communities. For the latter two, not only station coordinates and Earth orientation are of interest, but also the realization of the celestial reference frame (CRF), defined by positions of radio sources. VLBI has also been used to observe spacecrafts (e.g., Lebreton et al. 2005;Duev et al. 2012). Moreover, there are efforts to use the technique for the determination of planetary orbits (Jones et al. 2015) and relating those to the International Celestial Reference Frame (ICRF; Fey et al. 2015). Usually, differential VLBI ( VLBI) is used in order to cancel common error sources and achieve highly precise angular coordinates. Several experiments have been carried out with the purpose to determine the relative and absolute position of radio transmitters on the lunar surface as well as lunar satellites. The goal of these experiments was to obtain a better understanding of the physical properties of the Moon. Examples of recent lunar programs include LADEE (Lunar Atmosphere and Dust Environment Explorer; Elphic et al. 2014), CLEP (Chinese Lunar Exploration Program; Li et al. 2015) and GRAIL (Gravity Recovery and Interior Laboratory; Konopliv et al. 2013). The VLBI technique was also used for the spacecraft orbit determination during the SELENE (SELenological and ENgineering Explorer) mission (Kato et al. 2008) by using phase delays from same-beam differential VLBI (Kikuchi et al. 2009). The Chang'E-3 (CE-3) mission used the same-beam phase-referencing to obtain horizontal position estimates of the lunar rover w.r.t. the lander with a meter-level accuracy (Zhou et al. 2015). Another technique used in the course of this mission was the X-band two-way Doppler tracking, characterized by the measurement precision of about 0.2 mm/s (Wenlin et al. 2017). In terms of a surface position uncertainty, this can be translated 123 to an accuracy on the level of about 10 m (Williams et al. 2006).
The first VLBI observations of artificial radio sources on the Moon were carried out in 1970s using Apollo Lunar Surface Experiment Packages (ALSEP), deployed on the Moon during the Apollo program. Combination of S-band VLBI with range observations to laser retroreflectors allowed to determine the relative coordinates of the ALSEP transmitters, in the two transverse components, with an uncertainty of about 5 mas, which corresponds to about 10 m on the lunar surface (King et al. 1976). In the case of data obtained only through Lunar Laser Ranging (LLR), the retroreflectors' coordinates are known to a meter level, as stated by Williams et al. (2006).
Common to most of the aforementioned applications is that either astrometric VLBI or classical radio-astronomical mapping has been used to position the object of interest. Both approaches imply that positioning is only possible relative to a reference source and often a nearby phase calibrator is needed. In the latter case, relative measurements are also affected by position uncertainties of the reference source (Fey et al. 2015;Reid and Honma 2014). In order to avoid relative positioning or phase-referencing, one might try to position a target in an absolute sense, i.e., expressed directly in the celestial reference frame. This implies observing the source in geodetic VLBI mode and leads to the question how accurate and precise such an approach would be. In order to address this and related questions, we performed extensive simulations based on today's and next-generation geodetic VLBI schedules, to which we added observations of an artificial radio source on the Moon. In the following sections the design of the simulations and the obtained positioning performance is discussed and the implications on classical geodetic target parameters are studied.

Geodetic VLBI
Compared to Global Navigation Satellite Systems (GNSS), Satellite Laser Ranging (SLR) or Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS), geodetic VLBI is the only space-geodetic technique which allows to simultaneously determine all Earth Orientation Parameters (EOP; Sovers et al. 1998) and uniquely provides the Earth rotational phase (UT1-UTC). Geodetic VLBI is also important for the realization and the maintenance of the International Terrestrial Reference Frame (ITRF; Altamimi et al. 2016;Bachmann et al. 2016).
VLBI relies on pairs of telescopes which are separated by hundreds or thousands of kilometers and observe the same radio source simultaneously. The difference in signal arrival time (delay) and its time derivative (delay rate) are the main observables of geodetic VLBI and determined in a correlation process. By observing several different radio sources, one can estimate geodetic target parameters, i.e., station coordinates, EOP and source positions as well as auxiliary parameters (clock behavior, wet part of the troposphere; Haas and Schuh 1996;Behrend et al. 2000;Nothnagel et al. 2016). In today's geodetic VLBI, observations are carried out at two different frequency bands (X and S). This allows to compensate for ionospheric effects by forming a ionosphere-free linear combination. As the precision of the group delay measurements in geodetic VLBI is currently limited to several tens of picoseconds, there are efforts to improve it by developing the next-generation VLBI system called VLBI Global Observing System (VGOS; Petrachenko et al. 2009;Niell et al. 2006Niell et al. , 2014Cappallo 2014).

VLBI networks used in the study
In order to study the impact of observing a radio transmitter on the lunar surface along with classical geodetic VLBI sources, the so-called "rapid turnaround" IVS sessions were used for the following simulations. These sessions are named IVS-R1 and IVS-R4 and are observed on Mondays and Thursdays, respectively. The IVS-R1 sessions are scheduled based on a core network of 8-9 VLBI sites and, other than in the IVS-R4 network, the same core stations are used throughout the year. Thus, all IVS-R1 sessions of the year 2014 were chosen as the base for the following simulations in order to assess the average performance of the IVS-R1 network. A map of the VLBI stations participating in these sessions is presented in Fig. 1. In order to evaluate the concept of VLBI lunar observations in the upcoming VGOS era, a schedule from the system design phase was also used (Petrachenko et al. 2008;Behrend et al. 2009). Not only are VGOS observations expected to be more precise than present VLBI, but there will be also a significant increase in the total number of observations per session. Therefore, a network of 16 VGOS stations was considered ( Fig. 2) in order to investigate the benefits of combining lunar and quasar observations within future VLBI sessions. It has to be pointed out however that the VGOS network used in this study may differ from the one that is currently realized, in terms of quantity and distribution of stations. This means also that the twin telescopes as those in Norway (Spitsbergen), Sweden (Onsala) or Germany (Wettzell) were not considered in the following study. However, the existing VGOS stations occupy similar parts of the globe as those anticipated during the design phase of VGOS. In addition, the utilized number of VGOS-type observations per station is consistent with the current VGOS plans. Therefore, the reader should be able to identify, through this study, the benefits of lunar observations in the VGOS era, compared against the performance of the current VLBI systems (cf. Sect. 3.1).

Inclusion of lunar observations into geodetic VLBI schedules
For this study, lunar observations were added into geodetic VLBI schedules by modifying so-called VLBI experiment (VEX) files (Whitney et al. 2002). Those files are available for IVS-R1 sessions and can be downloaded via one of the IVS servers. 1 In the case of VGOS, we used a single VEX file converted from a SKD file 2 (Gipson 2010), which was created during the early design phase of VGOS. It originally consisted of 173,831 observations, which is about 40 times more than the average number of observations in a regular IVS-R1 session in 2014. In order to include lunar observations in the actual IVS-R1 schedules, we followed a very simple concept that did not 1 ftp://cddis.gsfc.nasa.gov/vlbi/ivsdata/aux/2014/. 2 ftp://gemini.gsfc.nasa.gov/pub/sked/net_sim/stat16_12_3p5D0ln. skd. require re-scheduling but simply replaced every n-th quasar scan with one to a radio transmitter on the Moon, regardless whether there was more than one scan at the same time. However, since such a target might not be visible from all stations that were scheduled in that particular scan, those stations at which the lunar target was not visible, i.e., below the local horizon, continued with their originally planned source. This concept was pursued for each IVS-R1 session as well as the chosen VGOS-type session. The study was carried out using different replacement strategies, by substituting every 2nd, 3rd, 4th, 6th and 8th scan in the original schedule (Table 1). Such modified schedules formed then the base for the actual Monte Carlo simulations, which allowed us to assess the impact of observing a transmitter on the Moon together with quasars. It has to be noted that there was a slight decrease in the total number of observations in schedules containing lunar observations, compared to the average number of 4,600 observations in the original IVS-R1 schedules. This decrease can be attributed to the replacement strategy used here where 123 The term replacement strategy refers to the frequency of substituted quasar scans within a single session some baselines are no longer part of a scan since only a few sites can see the Moon, whereas other follow the originally scheduled quasar. In the case of IVS-R1 schedules from the whole year 2014, the median number of baselines per scan was equal to three, both for lunar and quasar observations. For the VGOS schedule used in this study, the same statistics, for lunar observations as well as for all substitution strategies, amounted to ten. The incorporation of lunar observations for the VGOS schedule decreased the median number of quasar observations per scan from twenty-one to fifteen, only in the case of the substitution strategies S-2 and S-3. The fact that even a faint radio transmitter on the Moon will have a higher flux density than any of the natural radio sources considered in geodetic VLBI implies that the Signal-To-Noise-Ratio (SNR) targets ,which are implicitly set in each observing schedule, will not be violated when replacing a quasar with a lunar observation. This means that even if the slewing times would take longer to point the antenna toward the radio transmitter on the Moon, the scan itself would be shorter as the SNR targets are reached after a shorter integration time. Thus, it can be stated that this relatively simple replacement approach will lead to a feasible observation plan that allows to include a lunar target in a regular geodetic VLBI schedule. Moreover, in order to simulate a realistic location for a lunar-based radio transmitter, the source was assumed to be at 44.12 • N, 19.51 • W in the Moon's body-fixed reference frame. This position is close to the coordinates of a lunar lander and a rover which both were delivered to the surface of the Moon in late 2013 during the Chang'E-3 project .

VLBI delay models
Theoretical VLBI delays can be computed in accordance with the so-called Consensus model (Petit and Luzum 2010). When considering radio sources in the Solar System, this standard VLBI model is no longer valid and a more complex treatment for objects at a finite distance is required. Theoreti-cal models regarding near-field objects and the application of these models were described by, e.g., Moyer (2000), Klioner (2003), Sekido and Fukushima (2006) and Duev et al. (2012). In the latter model, which was used in the following study, the VLBI delay is obtained after solving the light-time equation at the reception time T 1 at the first station. The transmission time T 0 of the tracked spacecraft can be determined in an iterative way using T 1 (Duev et al. 2012;Moyer 2000): with R 01 referring to the distance between the barycentric position of the first station at its reception time T 1 and the position of the tracked object at the transmission time T 0 . A similar expression is used to solve for the reception time T 2 at the second observing station. The relativistic terms (RLT 01 , RLT 02 ) consider all planets of the Solar System as well as the Sun and the Moon. Hereby, the moment of the closest approach of the photon to the gravitating body or the position of the gravitating body at the retarded moment of time has to be included (Klioner 2003). Theoretical delays are obtained in the barycentric dynamical time (TDB) due to the fact that the motion of a spacecraft and VLBI telescopes are defined in the Barycentric Celestial Reference System (BCRS). Since the observed VLBI delays are referred the Terrestrial Time (TT), one needs to apply the Lorentz transformation in order to express the computed delay in the proper time scale (Sekido and Fukushima 2006). The total theoretical VLBI delay is then obtained by considering contributions from the atmosphere and technique-specific effects such as the antenna thermal deformation, antenna axis offset and displacement of the telescope reference point (Petit and Luzum 2010).
If one considers near-field VLBI delays for an artificial radio source located on the surface of the Moon, it is necessary to compute its position in the BCRS. This can be achieved by using planetary and lunar ephemerides such as the JPL ephemeris (Folkner et al. 2009). The latter contains also information on lunar librations, which are necessary to obtain the orientation of the principal axis (PA) system relative to the ICRS Earth equator and equinox (Archinal et al. 2011).

Simulation of quasar and lunar observations
All simulations were performed with the c5++ analysis software, which was developed with the goal to consistently process VLBI, GNSS and SLR (Hobiger and Otsubo 2014). This software is capable of generating simulated geodetic VLBI observations τ sim , which consist of a geometric delay τ g and contributions from the turbulent troposphere, reference clocks and thermal noise. In their most general form, such simulated observations can be expressed as where ε j denotes the elevation angle at station j. Values of zenith wet delay (ZWD j ) were simulated following the turbulence model of Treuhaft and Lanyi (1987), implemented in accordance with the algorithm described by Nilsson and Haas (2010) and mapped into slant directions with the wet mapping function (MF w ) introduced by Lagler et al. (2013). In order to simulate turbulent atmospheric conditions, stationdependent refractive structure constants and tropospheric heights were taken from Petrachenko et al. (2009) and wind speeds of 7 m/s in random direction, constant within a VLBI session, were assumed. The behavior of the VLBI station clocks (clk j ) was modeled as a sum of the random walk and integrated random walk processes (Herring et al. 1990) with a stability of 1 × 10 −14 for an averaging time of 50 min. The random observational error caused by the thermal noise (τ rnd ) was considered in the form of the Gaussian noise and added to the other contributions, as expressed in Eq. 2. The simulated atmosphere and the clock behavior correspond to realistic physical scenarios. The overall precision of the simulated delays is controlled by the standard deviation σ of the Gaussian random number generator, which produces values of τ rnd . The σ values are set station-dependent (e.g., σ 1 and σ 2 ) and error propagation according to σ = σ 2 1 + σ 2 2 is then used to obtain τ rnd . In the case of equally performing stations, the applied noise level per baseline is therefore √ 2σ j . Quasar observations were generated from a Gaussian random generator with standard deviation of 14.14 mm (47 ps), which corresponds to a station-specific thermal noise of 10 mm and reflects the average post-fit RMS of today's geodetic VLBI. For the simulations concerning the VGOS era (cf. Sect. 3.1), the baseline noise was set to 1.41 mm (4.7 ps).
The precision of VLBI group delay measurements depends on the SNR of the observation and the effective bandwidth obtained by synthesizing several few MHz wide channels, spread out over a much larger bandwidth. In the case of Xband signals emitted from the CE-3 lander, multiple tones separated by few tens of MHz are available for VLBI tracking. As stated by Zheng et al. (2014), this leads to an effective bandwidth of 38.4 MHz. The bandwidth synthesis technique was also utilized during the Chang'E-2 mission (CE-2) for precise orbit determination of the lunar orbiter. In accordance with Li et al. (2012), this yielded a precision of VLBI data on the order of 0.2 ns for an effective bandwidth of 20 MHz. Although, the signal strength of the CE-3 transmitter would imply very short integration times and group delay observables with a precision on the order of a few picoseconds, one needs to consider additional noise contributions and the fact that ionospheric delays remain uncalibrated (cf. Sect. 4). Hence, in order to study how target parameters are impacted by the precision of lunar observations, simulations for such observations were carried out with different Gaussian random noise levels. In total, 20 different levels for baseline noise were used in this study. They varied between 1.41 mm and 141.4 mm, spaced logarithmically.

Processing setup
In total, three different analysis options were tested. For each analysis option, 49 original and modified IVS-R1 sessions from 2014 were used. The modified schedules had lunar observations included following the replacement strategy described in Sect. 2.2. Each of these schedules was simulated 20 times for a certain level of lunar observation precision, expressed as corresponding Gaussian baseline noise. Moreover, each session was simulated 20 times based on the original schedule in order to derive a reference to which the performance of target and nuisance parameters could be compared to. A comparison of the reference solution against station position and EOP repeatabilities from the analyzed IVS-R1 sessions of 2014 confirms that the simulator performs close to real-world situations. This is also confirmed by the results presented by Kareinen et al. (2017), who made use of the same simulator for extensive studies concerning UT1-UTC determination.
Concerning the analysis, three different options were studied. Common to all options is the parametrization of the nuisance parameters ( Table 2). The first option, where station positions and EOP were assumed to be known and the two-dimensional (2D) lunar coordinates of the lander were estimated, was abbreviated L. When considering VGOS, this analysis option was named VGOS-L in order to indicate that only lunar lander coordinates were estimated, but using VGOS-type schedules and observations. The second analysis option, called LE, considered that the EOP were estimated as well. The third option (referred as LET ) was with the purpose to study a case where also all station positions are estimated together with the aforementioned parameters.
The coordinates of the lander were solved for as horizontal 2D components of lunar latitude and longitude (φ lan , λ lan ).  This parameterization was chosen since the VLBI observing geometry does not allow for a consistent and reliable decoupling of all three coordinates of the lander. Therefore, the height component was fixed to an arbitrary constant value of 0.00 m. However, results from additional simulations related to the three-dimensional case (3D), i.e., estimating all three coordinates of the lander, are presented in Sect. 3.2.
Earth rotation parameters, i.e., pole coordinates (x p , y p ) and UT1-UTC, were estimated session-wise. This applied also to station positions, in the case of the LET analysis option. Clock and troposphere parameters are estimated as described in Table 2. For all observations, an elevation cutoff angle of 3 • was applied. Since all station coordinates were solved for in the analysis option LET, the singularity of the normal-equation matrix had to be dealt with by introducing No-Net-Translation (NNT) and No-Net-Rotation (NNR) conditions w.r.t. the a priori ITRF coordinates (Petit and Luzum 2010).

Results
Repeatabilitites in the form of Weighted Root-Mean-Square (WRMS) errors were computed from the Monte Carlo results for each of the determined parameters. Since the WRMS refers to the a priori information, we can directly access the accuracy of each obtained parameter.

Position of the lander
We studied how the horizontal accuracy (WRMS 2D ) of the lander's position depends on the precision of lunar observa-tions. Figure 3 depicts this behavior for IVS-R1 sessions and substitution strategies S-2, S-3, S-4, S-6 and S-8.
As one would expect, more frequent and precise lunar observations lead to a higher accuracy of the lander's position. For observational precision worse than 20 mm, a direct relation between the noise of lunar observations and the WRMS 2D exists. On the contrary, observations with precision better than 20 mm do not improve the positioning results significantly and the WRMS 2D values settle between 270 and 340 mm. This is thought to be attributed to the fact that the tropospheric turbulence is dominating the noise budget and thus is the limiting factor, as discussed later in Sect. 3.5. The scatter plot as well as the corresponding histograms do not indicate any significant systematic effects.
Using single-frequency lunar observations, as in the case of the CE-3 probe, poses a problem for mitigating dispersive delays caused by the Earth's ionosphere. Assuming that ionospheric delay can be derived from Global Ionosphere Maps (GIM) (Schaer et al. 1996), for lunar X-band observations such corrections are expected to have an uncertainty of about 0.2 ns (Sekido et al. 2003;Hobiger et al. 2006). This corresponds to about 60 mm in terms of a baseline noise. Including ionospheric contributions into the error budget, one can now infer from Fig. 3 the expected accuracy of the lunar coordinates. This implies 2D horizontal accuracies around one meter.
In order to study how the estimation of additional parameters affects the lander's position, the WRMS 2D statistics were also computed for the LE and LET analysis options. We averaged WRMS 2D results from all five groups of modified schedules, i.e., from S-2 to S-8, which allowed us to compare the impact of different analysis options on the lander's In general, the obtained β factors indicate that the estimation of additional parameters decreases the accuracy of the estimated position of the lander. Depending on the precision of the lunar observations, this degradation is between 8 and 2% w.r.t. the L solution and indicates that highly precise lunar observations would be affected more than imprecise observations of the same type. Increasing the parameter space leads also to a larger uncertainty of the estimates. However, in the case of imprecise lunar observations this effect is negligible.
After studying the current VLBI performance in terms of lunar observations, the logical question is how the nextgeneration VLBI system will perform. VGOS will provide a measurement precision that is at least one order of magnitude better than the current geodetic VLBI technology. Thus, we followed a similar strategy as for the IVS-R1 sessions, but assumed an observation precision of 1.41 mm (4.7 ps) for quasar targets. The WRMS 2D values corresponding to VGOS-type schedules are presented in Fig. 5, including also a scatter plot of the lunar position estimates for the substitution strategy S-3 simulated with an assumed lunar observation precision of 15.97 mm. The pattern of the scatter is similar to the one shown in Fig. 3b and also indicates a larger spread in the latitudinal direction. The results show that VGOS-type observations improve the accuracy of the lunar coordinate estimates by a factor of eight compared to the IVS-R1 schedules (cf. Fig. 3). Although the impact of the troposphere is significantly reduced with VGOS, there is still a small contribution from this error source, thus preventing an improvement of the lunar lander's positions when decreasing the observational noise. This means that, for lunar observations with a precision of up to 10.0 mm, the WRMS 2D does not change significantly, varying between 20 and 50 mm in dependence of the number of lunar observations. For less precise lunar observations, the accuracy of the lander's position is still below a meter, which indicates the potential of lunar observations in the VGOS era. Given the poor sensitivity of geodetic VLBI observations w.r.t. all three position components of a radio transmitter on the surface of the Moon, only horizontal 2D positions have been considered here so far. Nevertheless, as shown in the following, in principle geodetic VLBI allows to determine coordinates of the lander in three dimensions. However, high correlations between position components decrease their precision and accuracy significantly. In order to study the capability of geodetic VLBI in terms of 3D positioning on the lunar surface, additional Monte Carlo simulations were carried out, in which also the height of the lander (h lan ) above the lunar reference surface was estimated. The L and VGOS-L analysis options were applied with this new parameterization.
The WRMS values of the estimated position components are shown in Table 3 for the lunar observation noise of 15.97 mm and based on the group S-3 of schedules.
When estimating also the height component of the lander, the position performance in terms of WRMS significantly decreases. In order to understand this behavior, one needs to consider the geometrical constraints that go along with Earth-based lunar observations. As discussed before, geodetic VLBI has a high sensitivity perpendicular to the line of sight and the direction of the baseline, but does not provide sensitivity in the direction toward a target on the Moon. As the Moon always faces the Earth with the same side, it is obvious that geodetic VLBI does not provide enough observations with suitable geometry allowing to decouple all three position components. This is reflected in the correlations (r ij ) between φ lan , λ lan and h lan derived from a single 24-h IVS-R1 session, see Table 4.
Since all three parameters are strongly correlated with |r ij | > 0.995, it is evident that it is almost impossible to solve for all three parameters unless external information is available or constraints are applied. On the contrary, when the height component is fixed, the correlation between the remaining two parameters φ lan and λ lan is significantly lower and amounts to − 0.685. Nevertheless, if all three position components need to be estimated from a single monolithic fit, the correlations prevent a meaningful interpretation of the individual coordinates. Individual components become indistinguishable and position solutions scatter within a spheroid that reflects the 3D accuracy of the lander's coordinates. In order to demonstrate this behavior, dedicated simulations with the IVS-R1 and VGOS-type sessions were carried out. As presented in Table 3, the W RM S 3D value for the VGOS is about three meters, which is about eight times smaller than the same measure obtained from IVS-R1 schedules. Even though the 3D positioning performance is worse than the 2D case for VGOS, as discussed in Sect. 3.1, the results are comparable to the vertical accuracy of present digital elevation models (DEM) of the Moon (Barker et al. 2016).

Earth orientation parameters
Having studied the potential of VLBI for positioning a radio transmitter on the surface of the Moon, one can now address the question how such observations impact classical geodetic VLBI target parameters. Among other estimates, EOP are products which are routinely derived from regular geodetic VLBI observations. Hereby, the LE and LET analysis options can act as proxies that give insights into this question. Therefore, the daily estimates of pole coordinates (x p , y p ) and the phase of the Earth's rotation (UT1-UTC) were studied based on the results from the Monte Carlo runs and compared against the performance of simulated IVS-R1 sessions which did not include observations to the Moon. Similar to the case of lunar coordinates, repeatabilities in the form of WRMS errors were used as evaluation criteria. Corresponding results are shown in Fig. 6. The results indicate slightly higher WRMS values for the x p and y p estimates when lunar observations are included in the IVS-R1 schedules. Also for UT1-UTC estimates a small accuracy degradation can be identified when including lunar observations. Although EOP estimates tend to be slightly worse than the results from reference solutions, the identified decrease in accuracy is rather small when comparing against the performance of actual VLBI sessions, especially in the case of UT1-UTC estimates. As stated by Malkin (2009), the precision of the latter parameter and pole coordinates, both derived based on the IVS-R1 sessions, amounts to 3 μs and 60 μas, respectively.
In general, the slightly worse performance can be explained by two main factors. First, due to the mutual visibility constraints, schedules with lunar observations contain less observations and thus lead to a slightly larger scatter of the EOP. Second, the inclusion of lunar observations, based on the replacement strategy presented before, leads also to less optimal observing geometries. Compared to the original IVS-R1 schedules, which have a varying observing geometry, lunar observations lead to a recurrent observing strategy due to a cyclic revisit of the lunar lander as a radio target. However, observations to artificial radio sources are characterized by short scan lengths. The latter feature paired with dedicated observing schedules could reduce the shortcomings of adding lunar observations to geodetic VLBI sessions. This indicates already how one could improve the concept in the future.

Station positions
With many quasar and lunar observations it is possible to examine to what extent the latter contribute to the determination of telescope positions. Hence, the LET analysis option allowed us to access the necessary statistics related to the station coordinates and their accuracies. Here, the WRMS values of position components expressed as topocentric coordinates (north, east, up) can be taken from the Monte Carlo simulations and understood as a measure for station position accuracy. As in the case of the EOP, WRMS values were computed for the original as well as the modified IVS-R1 observing schedules. Results for stations participating in more than 60% of the IVS-R1 sessions in 2014 and the substitution strategy S-3 are shown in Fig. 7. An insignificant increase in average WRMS of about 1 mm is found. Similar to the EOP-related study (cf. Sect. 3.3), this can be attributed to the smaller number of total observations per 24-h session, as compared to the original IVS-R1 schedules. In addition, one needs to take into account the change of the local sky coverage, which negatively impacts the sensitivity to estimate the station positions. These two drawbacks will be overcome when sophisticated scheduling strategies are considered.

Impact of error sources
Refractivity variations in the neutral atmosphere are considered to be the dominant error source which limits the accuracy of obtained VLBI target parameters and thus leads to a degradation of the estimated coordinates of the lander. This is especially pronounced in Figs. 3a and 5a where one can identify a settling of WRMS 2D for lunar observations with the precision of better than 10 mm.
The impact of the three simulated error sources (random noise, wet part of the troposphere, clock instability) on the accuracy of lunar 2D coordinates was investigated using the L option and S-3 schedules. This is summarized in Fig. 8. The first run considers all error sources and corresponds to the results shown in Fig. 3a. Assuming no contributions from the turbulent troposphere (while still solving for ZWD), as done in the second run, shows an improvement of about 50% in terms of WRMS 2D for lunar observation precision of better than 10.0 mm. Neglecting also the clock instability (while still solving for clk), leads to further improvement on the level of about 30%. In order to study almost noisefree observations, a dedicated simulation run assuming no noise contributions from the troposphere and clocks as well as VGOS-type quasar observations (i.e., σ = 1.41 mm) was carried out. This fourth run confirms that WRMS 2D linearly depends on the observation precision, as one would expect. In summary, it can be concluded that an improved handling of the tropospheric turbulence will lead to a better lunar positioning performance. Reducing contributions from the reference clocks and observational noise could lead to a centimeter position accuracy of the lander, as already partly confirmed by the VGOS study (cf. Sect. 3.1).

Summary and conclusions
We studied how geodetic VLBI can be used to determine the position of an artificial radio source on the Moon and we revealed the impact of this observational concept on classical geodetic parameters. We described the results of our extensive Monte Carlo simulations, based on the IVS-R1 schedules, and it was discussed how the quality and quantity of observations to the Moon affect the determination of EOP, VLBI station coordinates and position accuracy of the lander. We also carried out simulations reflecting the future VGOS performance. In addition, we investigated the capabilities to estimate three-dimensional lunar position components and revealed the limiting factors preventing centimeter-accurate positioning on the lunar surface.
In general, the results based on the IVS-R1 schedules (Sect. 3.1) indicate that a sub-meter position accuracy of an artificial radio source on the Moon could be achieved. Considering that ionospheric delays have to be corrected with the help of external information, such as GIM, a realistic positioning performance of about a meter is anticipated. However, if one would deploy a dual-frequency transmitter on the surface of the Moon, ionosphere-free observations would become available, leading to a horizontal position accuracy of better than 40 cm. Such performance could be already achieved when the lunar target is only observed in every eighth scan (i.e., 3% of total number of observations). Utilizing VGOS-type schedules provides an improvement in terms of the two-dimensional position precision by a factor of eight. This corresponds to WRMS 2D values in the range of 2-5 cm and 10-20 cm for ionosphere-free and externally ionosphere-corrected observations, respectively. Thus, for the VGOS case, the impact of the troposphere is no longer a major issue. Instead, the proper ionospheric delay mitigation for lunar observations is of high importance.
The results presented in Sect. 3.1 for the VGOS case are especially important from the scientific point of view. Positioning of an object in two dimensions with a precision of a few centimeters by geodetic VLBI can be beneficial for the development of lunar dynamical models and ephemerides. Geodetic VLBI could complement Lunar Laser Ranging (LLR) and contribute with sensitivity perpendicular to the line of sight. Moreover, long-term and frequent lunar observations could be provided within the IVS observing programs. Given that suitable lunar radio transmitters are provided, the existing global network of VLBI stations could track those as well as quasars within the same 24h experiments at no extra cost. As confirmed, inclusion of lunar observations into IVS-R1 schedules did not significantly impact the estimation of EOP and coordinates of the VLBI stations. Furthermore, observing programs, dedicated for both near-field and far-field observations, would allow to access the parameters that have been out of the reach for traditional geodetic VLBI, such as geopotential coefficients, lunar physical librations (Rambaux and Williams 2011) and parameters of general relativity theory.
VGOS-type quasar measurements have already been taken at several stations around the world. Moreover, the first observations of the CE-3 signals with geodetic VLBI were already performed in April 2014 on the Onsala-Wettzell baseline and regular observations of that probe with the global network of VLBI stations were also proposed to the IVS Program Committee which resulted in four OCEL sessions (Observing the Chang'E Lander with VLBI) per year from 2014 to 2016 Haas et al. 2017). This allowed to broaden the knowledge and gain new experience concerning geodetic VLBI observations of lunar targets. However, geodetic VLBI analysis of the Chang'E-3 data is still ongoing. In general, suitable processing routines for observations of different types of artificial radio sources (colocation satellites, GNSS; Tang et al. 2016;Plank et al. 2017) are still under investigation. Moreover, provision of geodetic VLBI observables suffers from technical difficulties which need to be solved before one can access the full potential of such new observation types .
The Chang'E-3 mission can be seen as a pathfinder for subsequent lunar missions with the aim to explore the Moon. A lunar lander equipped with space-geodetic instruments such as LLR reflectors and multi-frequency broadband radio transmitters would be a unique opportunity for comprehensive investigations concerning the structure and motion of the Moon and the dynamics of the Earth-Moon system (King et al. 1976). Geodetic VLBI can provide valuable insights that motivate the scheduling and observation of lunar radio transmitters along with standard quasar sources and therefore is expected to stimulate new observing concepts for geodesy as well as planetary science.