Baseline Vector Repeatability at the Sub-Millimeter Level Enabled by Radio Interferometer Phase Delays of Intra-Site Baselines

We report the results of position ties for short baselines at eight geodetic sites based on phase delays that are extracted from global geodetic very-long-baseline interferometry (VLBI) observations rather than dedicated short-baseline experiments. An analysis of phase delay observables at X band from two antennas at the Geodetic Observatory Wettzell, Germany, extracted from 107 global 24-hr VLBI sessions since 2019 yields weighted root-mean-square scatters about the mean baseline vector of 0.3, 0.3, and 0.8 mm in the east, north, and up directions, respectively. Position ties are also obtained for other short baselines between legacy antennas and nearby, newly built antennas. They are critical for maintaining a consistent continuation of the realization of the terrestrial reference frame, especially when including the new VGOS network. The phase delays of the baseline WETTZ13N – WETTZELL enable an investigation of sources of error at the sub-millimeter level. We found that a systematic variation of larger than 1 mm can be introduced to the Up estimates of this baseline vector when atmospheric delays were estimated. Although the sub-millimeter repeatability has been achieved for the baseline vector WETTZ13N – WETTZELL , we conclude that long term monitoring should be conducted for more short baselines to assess the instrumental effects, in particular the systematic differences between phase delays and group delays, and to find common solutions for reducing them. This will be an important step toward the goal of global geodesy at the 1 mm level.

measurements between pairs of nearby objects on the sky. Phase delays can be also used in geodesy to obtain relative positions between nearby antennas on the Earth with the highest accuracy, provided systematic effects are accounted for.
In the transition period of the geodetic VLBI systems, phase delays of short baselines enable significant scientific applications. Many antennas of the legacy VLBI system, which is mainly based on dual-band observations (2.3 and 8.4 GHz), though being continuously upgraded and still used, have reached the limits of their capability; this legacy system is pushed to the limits also because the Earth science studies continue to pursue more precise geodetic measurements. The next-generation geodetic VLBI system, known as the VLBI Global Observing System (VGOS; A. Niell et al., 2007;Petrachenko et al., 2009), has been developing worldwide with antennas of relatively small diameter, 12-13 m, and broadband receivers, 2.0-14.0 GHz, with the aim to achieve 1 mm station position accuracy and 0.1 mm/yr velocity stability on global scales. It is necessary to accurately tie these new, small antennas to the legacy, co-located antennas that have a long observing history since 1979 and have been playing a fundamental role in the realizations of the International Terrestrial Reference Frame (ITRF; Altamimi et al., 2016) to allow for a consistent continuation of the ITRF into the VGOS era. Recently, dedicated position tie measurements of this type have been performed, for instance, for the legacy antenna and the VGOS antenna at the Kokee Park Geophysical Observatory by A. E. Niell et al. (2021) and for the legacy antenna and the twin VGOS antennas at the Onsala Space Observatory by Varenius et al. (2021). An alternative way to derive these position ties is to make use of the global geodetic VLBI observations by the International VLBI Service for Geodesy and Astrometry (IVS; Nothnagel et al., 2017;Schuh & Behrend, 2012, see https://ivscc.gsfc.nasa.gov/index.html).
In this work, we analyze the observed phase delays to obtain position ties for as many co-located legacy and VGOS-compatible antennas and as many observations as possible. Our purpose is twofold: (a) to determine the baseline vectors between the legacy antennas and the co-located, new antennas with the highest possible precision and (b) to investigate the baseline vector repeatability of the short baselines determined from a time series of VLBI observations. The latter will allow us to separate the purely instrumental effects, affecting both short-baseline and long-baseline observables and dominating the estimates of the short-baseline vectors, from other contributions due to geophysical/astrophysical effects. The goals of this study are to contribute to the effort of the consistent continuation of the global Terrestrial Reference Frame (TRF) and to investigate the sources of error, mainly the instrumental effects, in VLBI observations.

Data and Data Analysis
We analyzed the IVS observations to derive the position ties for the antennas shown in Figure 1. The routine geodetic solutions of these global sessions have already been submitted by IVS analysis centers to the IVS combination center, which combines the results and provides the VLBI inputs for building the ITRF. (For the latest ITRF2020, the IVS analysis activities can be found at https://ivscc.gsfc.nasa.gov/IVS_AC/IVS-AC_ITRF2020. htm.) However, the short-baseline observables in these global geodetic VLBI observations can be analyzed independently from the observations of the entire network in each session in order to obtain the baseline vectors with the highest accuracy. The reasons are as follows: 1. In the routine geodetic VLBI solutions, observables at both S and X band are required to remove the dispersion affecting the radio signal when it passes through the charged medium, mainly the ionosphere. Any local radio interference, which is highly correlated for antennas at the same site, contributes large noise to the S band observables and thus to the ionospheric-free observables, though scaled down by a factor of 13.8. More importantly, false detections at S band lead to flagging the corresponding observables at X band as bad, and in not uncommon cases the observations of an antenna in one session are completely lost in the final data analysis due to the issues that happened only at S band. (See the comparison for baseline NYALES13S-NYALES20 in Section 3.1.2.) However, the observables at S band are not needed for short baselines, as the ionospheric effect is negligible for co-located antennas (pointing to a common source). 2. The position estimates of the co-located antennas treated independently in a geodetic solution of a full session are affected by systematic error sources, such as source structure, ionosphere, and atmosphere. In contrast, these systematic errors impose minimum impacts on the short-baseline observables. 3. Thermal noise can be one of the dominant errors in the short-baseline group delay observables, and it is significantly reduced by using phase delay observables (Ray & Corey, 1991).
4. Some of these short baselines are regularly scheduled in the VLBI sessions having a duration of 1 hr for the rapid determination of the highly variable Earth's rotation, called Intensive sessions, which by their design are not intended to be used for deriving station positions. They allow us to investigate the position precision that can be obtained from short-time observations, like the Intensive sessions.

Observations
In addition to the dedicated experiments for the Onsala antennas (the ONTIE sessions) that were reported in Varenius et al. (2021) and the Kokee antennas reported in A. E. Niell et al. (2021), short-baseline observations were found in three types of geodetic sessions: regular 24-hr sessions, special sessions of a combined network from legacy antennas and VGOS antennas, and Intensive sessions. The total number of the VLBI sessions (of these three types) analyzed in this study and the baseline lengths are reported in Table 1. Baseline ONSA13NE-ONSA13SW is formed by two VGOS antennas; each of the other 10 baselines consists of a legacy antenna and a new antenna with a small diameter in the 12-15 m range. The broadband receivers used in the VGOS system record the linearly polarized components of a signal, denoted by H and V, whereas the receivers of the legacy antennas are designed to record right-hand circular polarization, denoted by R. In the current data processing of VGOS observations, the pseudo-Stokes I visibilities are formed from the four linear polarization correlation products due to lack of knowledge of the cross-polarization "D" terms (see https://www. haystack.mit.edu/wp-content/uploads/2020/07/docs_hops_000_vgos-data-processing.pdf). For a mixed baseline of a legacy antenna and a VGOS antenna, a combined product of RH + RV is formed; the observations including these baselines in the network are referred to as the mixed-mode sessions. For the observations analyzed in this work, the first five baselines of Table 1 were observed in legacy S/X mode, and the remaining six baselines were observed in mixed mode. The new antennas involved in the first five baselines may have observed with a broadband receiver in other sessions or may be upgraded as a VGOS antenna in the future.

Phase Ambiguity
For geodetic VLBI observations in the legacy mode and the mixed mode, a multi-dimensional Fourier search from fringe phases of an interferometer gives (multi-band) group delay, delay rate, and visibility phase. The group delay is the derivative of phase with respect to frequency, whereas the phase delay is obtained as the ratio of the visibility phase to frequency. Phase delays Figure 1. Radio telescopes with position tie measurements reported in this study. At each very-long-baseline interferometry site (blue dot), there is a legacy telescope (black designator) and at least one new telescope (red designator). have intrinsically higher precision than group delays; however, they are typically not used in routine geodetic solutions due to unknown phase turns, that is, phase ambiguities.
Phase delay differs from group delay in terms of (a) the instrumental effects, such as the rotation of the feeds, the dispersion of the signal in the antenna system itself, and the signal delays in the waveguides prior to the injection of the phase calibration signals, (b) the frequency-dependent astronomical effects, the dispersive nature of the plasma along the line of sight and extended structure of radio sources, and (c) the magnitude of the thermal noise. The instrumental effects, which can be very large, either can be calibrated or are expected to be constant. The integrated plasma densities along the line of sight have very small differences for the co-located antennas of a short baseline. Most of the radio sources in the geodetic catalog, after a refinement over 40 years, are compact at the arcsecond scale, and the effects of structure for the majority of the sources at the scale of milli-arcsecond are relatively small for the short-baseline observations (see, e.g., Xu et al., 2016Xu et al., , 2019. For short baselines, the uncertainties of group delays due to thermal noise are generally far smaller than the phase ambiguity spacing. In the cases where the group delays are very noisy, for instance on the baseline NYALE13S-NYALES20, theoretical delays instead of the group delay observables can be used for directly aligning the phases over time, assuming that the unpredictable effects on the short baselines change relatively smoothly. Exceptional cases can happen when antennas have very unstable clocks or the a priori station position is very poorly known. The third option is to do a geodetic solution based on group delay observables for estimating the clock parameters and the station positions and then to employ them to connect the phases for resolving phase ambiguities.
The delay spacing of the phase ambiguities is about 120 ps at X band and about 450 ps at S band. In general, the variation of the differences between phase delays and group delays are expected to be relatively small compared to the ambiguity spacing, so that it is straightforward for most of the sessions to connect the phases over time.
Yet there can still be (ambiguous) constant offsets between phase delays and group delays, which will be fully absorbed by the estimated clock offsets. (We should note that resolving phase ambiguity is generally challenging for long baselines because of the impacts of, for instance, atmosphere.) In this study about short baselines, we used the group delays to eliminate the 2π phase ambiguity of the corresponding phase delays and afterward examined the differences between the group delays and the phase delays for all observations of a baseline in a session. If the differences over time follow the pattern of a smooth curve with a scatter significantly smaller than half the ambiguity spacing, it is an indication of successful elimination of phase ambiguities, while a failure would be obvious through a random distribution of the differences within the ambiguity spacing. This method was used as an initial inspection. The other methods were used as alternatives for some of the baselines.
When phase calibration signals are too weak to be useful (or not available) for removing the instrumental phase variations between various frequency channels, the observations of radio sources with high flux densities can be used as an alternative to calibrate the instrumental phases, which makes the fringe fitting of group delays possible. This process is referred to as manual phase calibration. However, in this case one may not be able to connect the phases because of the variations of instrumental phases over time. The details of the correlation process are written in the IVS correlator reports. The feed rotation angle (FRA) corrections need to be considered even for these very short baselines, since the two antennas at one site can have different mounting types leading to differences in the FRA corrections, as is the case for the two antennas HARTRAO and HART15M (equatorial/ altazimuth).

Comparison of Group Delays and Phase Delays
The differences between phase delays and group delays can be investigated after resolving phase ambiguities. These differences are shown in Figure 2 for four cases as examples, which demonstrate that phase ambiguities can be reliably resolved based on group delay observables.
There can be systematic variations in the differences, which can change as much as 100 ps over an hour, as shown for baseline WETTZ13N-WETTZELL in session 21MAY10XA. When estimating only a constant clock offset and a clock rate over the 24 hr (two parameters), the delay residuals from a solution of group delays in the session have a similar pattern as the differences between group delays and phase delays, whereas the delay residuals based on phase delays are much smaller and flat. The delay residuals are shown in Figure S1 in Supporting Information S1. This result strongly suggests that the differences are introduced by the group delays. They are largely absorbed by the clock parameters in a full geodetic solution. These effects may be caused by the dispersion effects in the waveguides of the receivers prior to the injection of phase calibration signals, the undesired wave reflections within the antennas, and spurious phase calibration signals (see, e.g., Rogers, 1991). Note that instrumental instabilities of this size will cause difficulties for resolving phase ambiguities of the observations on long baselines including one of these two antennas. This is one of the obstacles when resolving phase ambiguities for global geodetic VLBI observations and will be discussed in a future study. Such large variability occurs in other sessions including this baseline and in observations of other short baselines as well.
The systematic variations, though much smaller, are also visible on baseline HARTRAO-HART15M in South Africa and baseline ISHIOKA-TSUKUB32 in Japan. Recovering phase ambiguities for an Intensive session is shown in the bottom panel of Figure 2.
Based on closure analysis (see, e.g., Anderson & Xu, 2018;Xu et al., 2016Xu et al., , 2021, the inherently higher precision of phase delays can be seen directly at the observable level without a geodetic solution. Figure 3 shows the closure phase delays and closure group delays of triangle ONSA13NE-ONSA13SW-ONSALA60 in session ON0080 (20 March 2020). In principle, closures are sensitive only to the thermal noise and the effects of source structure, although the latter imposes minimum impacts on the observations of this small triangle for most of the geodetic sources. The unweighted and the weighted root-mean-square (rms) are 21.1 and 13.2 ps for the closure group delays, respectively, and they are 7.5 and 6.8 ps for the closure phase delays. Given that the thermal noise is independent among the three baselines, the noise level is about 12 ps in the group delays and 4 ps in the phase delays.
Considering that the dominating source of error in the short-baseline observables is the thermal noise, this improvement in the precision of observables can lead to significantly better results.

Ionospheric Corrections
The assumption that the ionospheric effect on short baseline is negligible can be validated after resolving phase ambiguities for both S band and X band observables. For baseline WETTZ13N-WETTZELL in session 21MAY10XA (a session with typical ionospheric delay corrections), the rms scatter of the ionospheric corrections at X band about the mean value, derived from the combination of the phase delays at S and X band, is less than 1 ps. For baseline NYALE13S-NYALES20 (in session 21JUN24XE), about 1.5 km apart, the rms scatter is 2 ps. For baseline ISHIOKA-TSUKUB32 (in session 16DEC20XA), about 16.6 km apart, the rms scatter is similar to the values for the baseline NYALE13S-NYALES20.
In order to assess how much the S band observables corrupt the short-baseline observables in routine geodetic solutions, the rms scatters of the ionospheric corrections at X band, derived from the group delays at S and X band in the conventional way and restored in the databases, are calculated for the short baselines in the mixed mode session RD2005. The rms scatter is 15 ps for baseline WETTZ13S-WETTZELL of 0.2 km length and is 90 ps for baseline ONSALA60-ONSA13SW of 0.5 km length. With this justification, the phase delays at S band were not used in our solutions because they can lead to flagging as outliers a significant amount of useable X band phase delays.

Data Analysis
In the multiple steps of VLBI data processing (see https://ivscc.gsfc.nasa.gov/about/resolutions/ IVS-Res-2019-02-AnalysisLevels.pdf), geodetic analysis is performed with the aim of estimating the parameters of geodetic interests, such as Earth orientation parameters (EOP), station positions, and source positions. In the geodetic solutions of short-baseline observations, there are two other possible kinds of parameters in addition to baseline vectors: clocks, accounting for the relative behaviors of the two frequency standards and for the instrumental delays, and differential zenith wet delays (dZWDs), accounting for the different atmospheric effects between the two antennas. The clocks were characterized by a continuous piece-wise linear (PWL) function with a time interval of usually one hour (but half an hour for baseline SESHAN13-SESHAN25, see Section S1.1 in Supporting Information S1). For the atmospheric delays, the hydrostatic part was modeled, while the impact of the wet path delays due to water vapor was investigated by comparing the results of the baseline vector estimates from not estimating dZWDs Figure 3. Comparison of closure group delays (blue squares) and closure phase delays (red squares) for the triangle ONSA13NE-ONSA13SW-ONSALA60 in session ON0080 demonstrates the significantly higher precision of the phase delays than the group delays. The scatter of the closures indicates the contributions of thermal noise. The red squares marked by black circles indicate the observations of source 3C274, which is well known to have large scale structure and has a similar pattern in its closure phases from the other ONTIE sessions with these three antennas. The closures suggest that about 2% of the short-baseline observations may be significantly affected by source structure at large angular scales. and estimating them using PWL functions of different time intervals (see Section 5 (9) for discussion). The results reported in Tables 2 and 3 are obtained for the two baselines ISHIOKA-TSUKUB32 and NYALE13S-NYALES20 (longer than 1.0 km) by estimating dZWDs with a continuous 1-hr interval PWL function and for the remaining baselines by not estimating dZWDs. Geodetic analysis was carried out by using either phase delays or group delays.
The software package νSolve (open source, available at https://sourceforge.net/projects/nusolve/) was used for the geodetic analyses. For each session, we reset the configuration in the original databases to remove the flagging and weighting information and the ionospheric corrections, excluded the observations of all antennas apart from the two antennas of the desired short baseline in a session, restored all the useable observables, examined and adjusted the phase ambiguities in a program developed by ourselves, flagged the outliers, and performed the solution based on either group delays or phase delays at X band. In geodetic solutions, as guided by the νSolve user manual, one step that is commonly used is to determine a baseline-dependent uncertainty in addition to the formal error of each observable in order to derive a more realistic error used for the weighting; this additive uncertainty σ add is a constant value in a session for each baseline and is determined in an iterative way until the reduced χ 2 is unity.

Results
When there are more than three sessions available for a baseline, the baseline vector repeatability is defined as the weighted root-mean-square (WRMS) scatter of the relative position estimates from these multiple sessions about the weighted mean value. We evaluated this metric for the three components of a baseline vector and present the results always in the sequence of the east, north, and up directions.  Note. The topographic coordinate system in this work is defined to be centered at the position of the first antenna of a baseline. The baseline vector WETTZ13N-WETTZELL from local survey is reported for comparison. L is baseline length with uncertainty σ L . a Baseline length L is derived as the mean of the baseline length estimates over multiple sessions in the same way as for the three position components; therefore, there can be a discrepancy of a few tenths of millimeter between the reported L and the value that one can calculate from the root of the sum of the squares of the three position components. Uncertainty σ L is calculated from the time series of the baseline length estimates as the uncertainty of the mean value instead of doing error propagation from the uncertainties of the three position components. This process provides an evaluation of baseline length as an independent quantity and is used in the study for these four baselines.

WETTZ13N-WETTZELL
Geodetic/astrometric VLBI makes routine observations of tens of radio sources typically for 24 hr or for 1 hr in one session. These two antennas have simultaneously participated in these two types of IVS observations since 2015 (Schüler et al., 2015).
The correlator centers for processing VLBI observations by using the fringe fitting program fourfit started to apply a special mask called notch filter to mitigate the corruption due to specific phase calibration signals after October 2018. The width of such a notch filter depends on the spectral resolution which is used for correlation: the higher the resolution, the narrower the notch filters. Therefore, only the sessions since R4889 (11 April 2019) or processed by the correlators after 1 May 2019 have useable observables on baseline WETTZ13N-WETTZELL. There are 107 sessions as listed in Table  S1 in Supporting Information S1. (We note that reprocessing the observations of this baseline since 2015 from visibility data will produce four more years of useable observations.) The results from the analyses of not estimating dZWDs for this baseline are presented here, whereas the results of estimating them will be discussed in Section 5. The mean number of total observations in these 107 sessions is 302, and the mean number of used observables in the solutions is 276 and 277 for group delay analyses and phase delay analyses, respectively. The mean value of the WRMS delay residuals is 15.6 ps for group delay analyses and 3.9 ps for phase delay analyses. They are approximately at the same level as those determined by the closures of the triangle formed by the Onsala antennas.
The mean formal errors of the estimates of the baseline vector in the east, north, and up directions are 0.6, 0.6, and 1.3 mm from group delay observables, respectively, and they are 0.2, 0.1, and 0.3 mm from phase delay observables. Because the estimates from different sessions are scattered more than one would expect from their formal errors, the formal errors were inflated by introducing a constant additive uncertainty such that the reduced χ 2 of the time series of each coordinate component becomes unity. The additive uncertainty is an indication of the systematic error level in the results that is not measured by the (original) formal error. They are 0.8, 2.5, and 2.3 mm for the three position components from group delay analyses, and 0.3, 0.3, and 0.7 mm for phase delay analyses, respectively. The results suggest that the sub-millimeter precision can be achieved for all the three components of this short baseline by phase delays in a single 24-hr session with the S/X observing mode.
We used the inverse of the sum of the squares of the formal error and the additive uncertainty as the relative weight for each individual estimate from one session in calculating the weighted mean baseline vector and the repeatability. The weighted mean of the baseline vector estimates from both group delays and phase delays are presented in Table 2. The baseline vector repeatabilities are 0.3, 0.3, and 0.8 mm from phase delay analyses and 1.1, 2.6, and 3.0 mm from group delay analyses. The precision obtained for this 123 m baseline based on phase delays is likely to demonstrate the best performance that the geodetic VLBI system with the S/X observing mode is capable of. The repeatabilities of the position of WETTZELL based on group delays in the 24-hr global sessions are 3, 5, and 9 mm according to the IVS internal report of the ITRF2020 on the 20th IVS analysis workshop in 13 September 2021.
The residuals of the estimated Up coordinates of antenna WETTZ13N from both phase delays and group delays are shown in Figure 4. There is a significant difference in the up direction between group delay and phase delay results; the weighted mean of the differences in the Up component (in the sense of group delay result minus phase delay result) is −1.7 ± 0.2 mm. The distribution of the residuals in the horizontal plane is shown in Figure 5. The majority of the east and north residuals from the phase delay analyses are within ±0.5 mm. The residuals from group delay analyses systematically spread in the north direction, but they do not show a temporal dependence. The results from group delays in the 24-hr global sessions also have a larger scatter in the north direction than in the east direction as shown in the IVS internal report of the ITRF2020. The differences in the mean horizontal components between group delay results and phase delay results are within three times the uncertainties of the group delay results. Complementary to the 24-hr sessions, Intensive sessions have been carried out since 1984 (Robertson et al., 1985) to rapidly determine Earth's highly variable phase of rotation. They last for one hour and currently are observed every day by two globally spaced antennas, generally WETTZELL and KOKEE, and every Monday by more than two antennas including WETTZELL and WETTZ13N. Due to continuous improvements in VLBI antenna sensitivities and in scheduling (see, e.g., Baver & Gipson, 2020;Schartner et al., 2021), and taking advantage of the consequent more even distribution of useable radio sources on the sky, it has become possible to estimate relative positions for the co-located pair WETTZELL and WETTZ13N from the Intensive sessions. We use the Intensive sessions here to learn how well the short baseline vector can be determined from one-hour observations by comparing to the results obtained from 24-hr observations.
We have processed 58 Intensive sessions that included the WETTZ13N-WETTZELL baseline within the same time period as the 24-hr sessions, listed in Table S2 in Supporting Information S1. These Intensive sessions on average consist of 43 scans, and the mean number of useable observables of the baseline is 39.0 and 39.6 for group delay and phase delay, respectively. This is at least twice the average number per hour of 24-hr sessions. In the data analysis of the Intensive only six parameters are set up: three for the baseline vector and three for the clock. The means of the WRMS delay residuals are 9.5 ps for group delay analyses and 3.1 ps for phase delay analyses, which are significantly smaller than those of the solutions based on the 24-hr sessions. However, since the Intensive sessions do not have significantly higher signal to noise ratio than the regular IVS global 24-hr sessions, these smaller delay residuals in the Intensive sessions may indicate that the shorter sessions are over-parameterized.
The residuals of the position estimates from the Intensives with respect to the reference position obtained from 24-hr sessions are shown in Figure 6. The means of the residuals from phase delays are −0.19 ± 0.16, −0.23 ± 0.14, and 0.20 ± 0.17 mm in the east, north, and up directions, respectively; they are 0.07 ± 0.44, −0.19 ± 0.51, and 0.71 ± 0.60 mm for group delays. The formal errors of the position estimates based on phase delays are on the level of 0.6, 0.6, and 0.8 mm for the three components, respectively; the additive uncertainties to these formal errors are 1.0, 0.8, and 1.0 mm in order to get the χ 2 of the residual time series being unity. The differences in the mean positions between the 1-hr observations and the 24-hr observations are not significant with respect to their uncertainties. However, the position residuals show systematic variations, mainly in southwest and northeast as shown in Figure 6. The phase delay analyses produce a baseline vector repeatability of 1.3, 1.1, and 1.3 mm, and the group delay analyses result in 3.4, 4.0, and 4.5 mm. Phase delays on a short baseline in an Intensive session have a capability of determining baseline vectors at the 1 mm level.

NYALE13S-NYALES20
The legacy antenna NYALES20 in Norway has an observing history of about 30 years, and it is still one of the most active geodetic stations. The new antenna NYALE13S has participated in the IVS sessions since early 2020 and operated through a series of shakedown experiments; the legacy antenna NYALES20 observed many sessions in 2020 and 2021 with a warm receiver. Thus, the observations of this baseline often have large noise contributions. Due to the large measurement noise and the poor a priori position of the new antenna, it can be challenging to eliminate the phase ambiguities for this baseline. We have 19 sessions available for this baseline to perform both phase delay and group delay analyses. The atmospheric effects were modeled as a PWL function with an interval of one hour, because the impact of wet delays is expected to be larger for this 1.5 km baseline than for baseline WETTZ13N-WETTZELL and the current data do not allow us to have a clear conclusion for the case of NYALE13S-NYALES20. The number of used observables and the WRMS delay residuals based on the two types of observables are reported in Table S3 in Supporting Information S1. The mean of the WRMS delay residuals from the IVS reports of the routine data analysis, labeled as "S/X band delays" in the table, is 53 ps for on average 220 used observables, and the residuals are significantly larger than the typical measurement noise level in the geodetic observations. By removing the involvement of the S band observables in the analyses, labeled as "Group delays" in the table, the number of useable observables increased by 34%, and the mean of the WRMS delay residuals decreased to 35 ps. A significant improvement has been obtained by using group delays at X band. The mean of the WRMS phase delay residuals is about 16 ps, a significant reduction from the group delay value.
The weighted mean estimate of the baseline vector calculated as the mean estimate is reported in Table 2. The phase delay results yield a baseline vector repeatability of 1.3, 1.1, and 4.7 mm, and the group delay analyses yield 2.0, 1.8, and 16.0 mm. The repeatability determined by group delays in the up direction is one order of magnitude larger than that in the horizontal directions, which means that either the large noise level in group delays has a larger impact on the up direction or the large noise in the group delays are not purely random but includes some systematic errors. As the group delays of this baseline obtained from manual phase calibration have a significantly lower noise level than those based on the observed phase calibration signals (as seen in the session 21JUN17XE, which was reported both ways), the issues in the group delay results may be due to the phase calibration systems. Referred to the mean position from phase delays, the residuals are shown in Figure 7 for both phase delays and group delays. The difference between the mean estimates from group delay analyses and phase delay analyses is within the uncertainties of the group delay results in the horizontal plane and about four times the uncertainty in the up direction.

ISHIOKA-TSUKUB32
The new antenna ISHIOKA can observe with both the S/X and the VGOS mode, thanks to the interchangeable S/X band and broadband receivers. The 32 m antenna TSUKUB32 ended observing in December 2016, and it was dismantled in 2017. It is only possible to derive the position tie for this 16.6 km baseline by analyzing historical observations. The reported results of this baseline are from estimating the dZWDs with an interval of 1 hr. The phase delays in the 17 sessions listed in Table S4 in Supporting Information S1 produce a baseline vector repeatability of 0.9, 0.9, and 3.4 mm, and the group delays give 1.2, 2.0, and 4.7 mm. The reference position and the mean position from group delay analyses are listed in Table 2. The difference between the group delay and phase delay results is 8.2 mm in the up direction, significant at the 5-sigma level, and is about 1 mm in the horizontal plane.

HARTRAO-HART15M
Because the frequency standard of antenna HART15M was tuned down by ∼4.5 Hz, this short baseline has useable observables without applying the notch filter in correlation. The solutions based on the phase delays of this baseline are only slightly better than the ones based on group delays as indicated by the WRMS delay residuals in Table S5 in Supporting Information S1. The mean estimates of the baseline vector from phase delays and group delays in the eight sessions are presented in Table 2. The differences between the results from the two types of observables are not significant with respect to their uncertainties.  Table 2. Note the different scale compared to the residual scatter shown in Figure 5.

Baselines With Only One or Two Global Sessions
The results of the baseline vectors from phase delays for the baselines with only one or two global sessions available are reported in Table 3. The VLBI data for the baseline SESHAN13-SESHAN25 are from session AOV056 (3 February 2021), and the data for the other six baselines are from two mixed sessions, RD2005 (24 June 2020) and RD2006 (8 July 2020). The detail of the data analysis of these observations are presented in the Supporting Information S2.

Comparison of the Results
Local survey measurements have been carried out at the Wettzell site to obtain the baseline vectors with an uncertainty of about 0.5 mm in each of the three components (Mähler et al., 2019). The local-survey result of the baseline vector WETTZ13N-WETTZELL is reported in Table 2. This result was derived from measurements over the course of several years, and thus has no nominal temperature of the local environment. Nevertheless, as we will see in Section 5, the baseline vectors among the antennas at the Wettzell site are very insensitive to the thermal expansion on the three antennas. The difference of the baseline vector from phase delays with respect to the local survey is −0.5, −0.3, and −1.6 mm in the east, north, and up directions, respectively, and it is −0.4, 0.3, and −3.3 mm for group delays. The VLBI results do not significantly differ from the local-survey tie in the horizontal directions, but the differences in the up direction are significant. The Up component from phase delays is closer to the local survey than from group delays for this baseline. As reported in Table 3, the local-survey tie of another short baseline at Wettzell, WETTZ13S-WETTZELL, has a significant difference (at the 3-σ level) in the east direction with respect to the VLBI results, about 1.5 mm. The comparisons of these two sessions suggest that the VLBI results and the local survey measurements may have a difference of 1-2 mm in the horizontal directions and up to a few mm in the up direction. However, we cannot draw a similar conclusion about the results of other baselines due to unmodeled systematic errors, as will be discussed in the next section.
As mentioned, the results of the short baselines at the Kokee Park and at the Onsala Space Observatory were previously reported by A. E. Niell et al. (2021) and Varenius et al. (2021), respectively. Compared with the result of the baseline KOKEE12M-KOKEE from VLBI measurements with a mean date of 11 April 2016 (A. E. Niell et al., 2021), the difference between our results is insignificant in the horizontal directions but has a magnitude of 4.0 mm in the up direction, which may be due to the lack of phase cal correction for the KOKEE antenna, as noted in A. E. Niell et al. (2021). For the legacy antenna and the twin VGOS antennas at the Onsala Space Observatory, the difference in the baseline vector ONSA13NE-ONSA13SW between our results and that reported in Varenius et al. (2021) is less than 0.2 mm in the horizontal directions and 0.8 mm in the up direction.
As an independent evaluation of the baseline vectors, our results from phase delays were compared to the latest realization of ITRF, ITRF2020 (see https://itrf.ign.fr/en/solutions/ITRF2020). See Table 4. Two baselines not listed in the table: ISHIOKA-TSUKUB32 because of the post-seismic deformation model employed for station TSUKUB32 in the ITRF2020 but not for station ISHIOKA, and SESHAN13-SESHAN25 due to the missing station SESHAN13 in the ITRF2020. The nine baselines listed all have position differences larger than 1 mm, even though most of them are consistent within the uncertainties that are dominated by those of the ITRF2020. Two baselines have differences at the cm level: NYALE13S-NYALES20 and RAEGYEB-YEBES40M. The former baseline vector has an uncertainty of several centimeters in the ITRF2020 due to the large impact from S band as discussed in Sections 2.4 and 3.1.2, and the latter one is most likely affected by the receiver replacement at YEBES40M in 2011.  Note. L is baseline length with uncertainty σ L . The results are from geodetic solutions using phase delays. The baseline vector WETTZ13S-WETTZELL determined from the local survey measurements is reported.

Discussion of Sources of Error
It is worthwhile to note that sub-millimeter repeatability has been demonstrated through short baselines by other space geodetic techniques than VLBI, for instance, the Global Navigation Satellite Systems (GNSS) by Hill et al. (2009); King and Williams (2009). These studies have used the short-baseline time series to investigate the site-specific errors and the stability for GNSS. However, due to completely different data collection and processing methods between the various space geodetic techniques, the sub-millimeter repeatability and the error investigation should be carried out independently for each technique toward the 1 mm accuracy goal of global geodesy.
An early phase delay analysis of 11 VLBI sessions of the 1 km baseline HAYSTACK-WESTFORD determined a baseline vector repeatability of 5, 3, and 7 mm in the east, north, and Up components (Carter et al., 1980). Then, later with the improved Mark III VLBI system, Herring (1992) obtained a repeatability of 0.8, 0.7, and 2.3 mm for the same baseline by using phase delays in 24 VLBI sessions. The baseline length of WETTZ13N-WETTZELL is 123 m, one order of magnitude shorter than that of baseline HAYSTACK-WESTFORD; both baselines observed in the S/X mode. As a continuation of the investigations from six sessions of WETTZ13N-WETTZELL in Halsig et al. (2019), which were based on group delays, the larger data set of 107 global 24-hr sessions and 58 Intensives over 2.5 years and the significantly improved repeatability in our study provide the opportunity to assess the error components more stringently.
What are the important sources of error in the repeatability of the baseline vector estimates of this short baseline?
The investigation based on the metric of repeatability in most cases is sufficient for geophysical and astrophysical studies, such as a change in the orientation of the Earth in space and the station position variations due to tidal displacements or tectonic motions. However, some of the instrumental effects can be highly repeatable and thus are not detected by the WRMS scatter of the estimates. It is necessary to also investigate those repeatable errors in the VLBI system itself for the purpose of combining the station positions from various space techniques for the realization of ITRF. Therefore, which of these error components are repeatable (related to accuracy) or not repeatable (reliability and precision)? We devote this section to addressing these questions.
(1) Measurement noise in delay observables. The uncertainty of a baseline vector estimate due to measurement noise in observables generally manifests itself as the formal error from the geodetic solution based on least square fitting (LSF). A session-wise delay noise was added in quadrature to the uncertainties of the observables in LSF to account for potential systematic errors already at the observation level. And in fact, the additive noise is comparable to the corresponding WRMS delay residual as shown in Table S2 in Supporting Information S1. Therefore, the uncertainty due to measurement noise is lower than the formal error of the estimate by a factor of approximately 1/ √ 2 . Taking this factor into account, the impact of measurement noise in phase delays on the baseline vector of WETTZ13N-WETTZELL is 0.1 mm on the horizontal plane and 0.4 mm in the up direction, and for group delays it is 0.5 and 2.0 mm.
(2) Cable delays. The time delays of the astronomical signal passing through the electronic devices within the antennas are expected to be smoothly varying and are absorbed in the clock parameters. In order to eliminate their variations, phase calibration systems are used to correct the visibility phases. Meanwhile, the time delays through the cable that carries the precisely timed pulses from the frequency standard to the injection of phase calibration signals are actually not experienced by the astronomical signals; these cable delays are measured as corrections to be applied in geodetic solutions. Of the 11 baselines in this study, only KOKEE12M-KOKEE, RAEGYEB-YEBES40M, WETTZ13S-WETTZELL, and the baselines formed by the Onsala antennas have the cable delay corrections available for both antennas. In fact, proxy corrections for KOKEE12M instead of direct measurements were used (see, A. E. Niell et al., 2021), and the cable delay corrections of antenna RAEGYEB were not applied for the final solution. It is possible that variations of the time delays in the antenna electronics and cabling cause significant impacts due to the missing corrections for some antennas.
We ran geodetic solutions of the 107 sessions in which the cable delay corrections of antenna WETTZELL were not applied by intention. We must emphasize that turning off the cable delay corrections did not degrade the solutions in the sense of the WRMS delay residuals and the baseline vector repeatability when comparing to the solutions with the corrections applied. After turning off the cable delay corrections of antenna WETTZELL, the mean of the changes in antenna WETTZ13N Up estimates from phase delays is −0.23 ± 0.03, −0.07 ± 0.02, and −1.65 ± 0.11 mm in the east, north, and up directions, respectively. As discussed in Section S1.3 in Support-ing Information S1, the cable delay corrections of antenna RAEGYEB have an impact of about −5 mm on the Up component of RAEGYEB if the cable delay corrections are not applied. For the ONSALA60 antenna, the impact is estimated to be dominant also in the up direction, which is even at the level of 1 cm (Varenius et al., 2021). The uncertainty due to missing cable delay calibrations may be at the sub-millimeter level on the horizontal directions and a few millimeters or larger in the up direction. The impact is repeatable at the 0.1 mm level for antenna WETTZELL. This repeatable feature may affect other antennas if the distribution of the elevation and azimuth angles does not change dramatically from session to session such that the cable is twisted in the same manner.
(3) Antenna thermal expansion. The antenna structure experiences thermal deformation, and this leads to station position changes due to the temperature variations at the site. This effect and the models have been well studied (see, e.g., Nothnagel, 2009, and the references therein). The antenna-dependent parameters in these models are maintained and publicly available for most of the geodetic antennas (from https://raw.githubusercontent.com/ anothnagel/antenna-info/master/antenna-info.txt). However, the practical problem of applying these models at the observation level in geodetic solutions can be that the desired temperatures of the antenna structural elements are specific functions of the time history of the ambient temperatures. Therefore, the information of these temperatures are not complete in VLBI databases. Nevertheless, the mean temperature during the 24 hr of observations would not significantly differ from the mean value of the temperatures that actually cause the thermal expansion. We use the temperatures recorded at the same epochs of observations in a session to derive the mean value and assess this source of error accordingly.
The temperature dependence of the contribution of thermal expansion to the raypath length consists of constant terms and elevation-dependent terms. In the case of WETTZ13N and WETTZELL, which have zero axis offsets, the elevation-dependent terms are proportional to the sine of the elevation angle, and the extra path length due to these terms exactly mimics the effect of a vertical displacement of the antenna. The antenna dimensions of the elevation-dependent terms are: (a) −2 m for the relative foundation height of WETTZ13N minus WETTZELL with an expansion coefficient of 1.0 × 10 −5 per °C and (b) +1 m for the relative height of the supporting axes with an expansion coefficient of 1.2 × 10 −5 per °C. Therefore, an increase of 1°C at the site will cause the Up component of WETTZ13N's position (relative to WETTZELL) to decrease by 0.008 mm. Based on the meteorological data recorded at each observation epoch, the median value of the mean temperatures in these sessions is 7.7°C, which is very close to the reference temperature of the model parameters previously mentioned. The rms scatter of the temperatures within one session has a median value of 2.3°C; thus, the impact due to the temperature variations within a day is negligible. There is a seasonal variation in the mean temperatures with an amplitude of 11.5°C ± 0.7°C. Even though the actual variations due to this effect will be about −1.4 mm in winter and about +1.4 mm in summer for both antennas, the corresponding impact on the baseline vector of WETTZ13N-WETTZELL has a magnitude of only 0.1 mm in the Up component.
For baseline ONSALA60-ONSA13SW from 24 June to 8 July 2020, the impact of thermal expansion is about 0.2 mm in the up direction due to the temperature change of 5.2°C at the site. We also note that the temperatures can be different between these two antennas because ONSALA60 is enclosed in a radome. Nevertheless, the impact of this effect cannot explain the difference in the estimates of the Up component of this baseline vector from the two sessions, as presented in Table 3.
The impact of thermal expansion is typically very small (i.e., below 1 mm) for the short baselines except RAEGYEB-YEBES40M, since the temperature is close to the same for multiple antennas at a site and since they tend to have similar physical dimensions. This source of error should show a seasonal variation in the up direction.
(4) Antenna gravitational deformation. In contrast to the constant, known coefficients of thermal expansion and antenna structure, the gravitational deformation of the main reflector and its distance to the secondary reflector need to be measured individually for each antenna. Extensive measurements of this effect for the Onsala antennas have been made by employing other measuring techniques (see, e.g., Bergstrand et al., 2019;Lösler et al., 2019;Nothnagel et al., 2019). These studies demonstrated that the effect produces systematic offsets in the Up component of the position of the 20 m Onsala antenna of about 1 cm; however, the change is smaller than 1 mm for the 13.2 m VGOS antennas at the site. These results may be considered representative of the impacts on the legacy antennas and the VGOS antennas at the other geodetic VLBI sites. The IVS recently adopted the resolution of every radio telescope operating in IVS observing sessions being surveyed for gravitational defor-mation investigations (see https://ivscc.gsfc.nasa.gov/about/resolutions/IVS-Res-2019-01-TelescopeSurveys. pdf).
The extra raypath introduced by antenna gravitational deformation is a smooth curve with respect to elevation angle, as measured for ONSALA60, similar to a sine function. Therefore, the major impact on station position estimates is on the Up component, as shown in the data analysis of Varenius et al. (2021). If this elevation dependence were exactly a sine function, the impact would again be equivalent to a displacement in the Up position with a magnitude of the gravitational deformation at zenith direction. Moreover, regardless of the exact form of the elevation dependence, the impact on position estimates is repeatable as long as the elevation angles have a similar distribution from one session to another. This is valid for baseline WETTZ13N-WETTZELL, because there is no significant difference in the scheduling of the 107 sessions with the same goal-only four of the 107 sessions are for the determination of TRF and the other 103 are R1 and R4 sessions. Nevertheless, we might conclude that the non-repeatable component of the antenna gravitational deformation is likely smaller than the repeatability of the baseline vector WETTZ13N-WETTZELL and thus the omission of gravitational deformation correction doesn't affect our assessment of repeatability. It would affect the accuracy of the results.
Since the gravitational deformation of the VGOS antennas is measured to be below 1 mm (though to be measured and confirmed for more VGOS antennas), an order of magnitude smaller than that of the legacy antennas, the corresponding error for the VGOS network is expected to be smaller than that of baseline WETTZ13N-WETTZELL.
(5) Antenna tilt. It is possible that the supporting axis of an antenna tilts toward a direction in the horizontal plane as time goes on (A. E. Niell et al., 2021). Based on the 3.5 years of observations of WETTZ13N-WETTZELL, however, the horizontal motion is determined to be negligible, −0.02 ± 0.04 and −0.03 ± 0.03 mm/yr in the east and north directions, respectively.
(6) Signal chain. Due to the large systematic differences in group delay and phase delay observables as shown in Figure 2, we believe that there is a possibility of significant errors introduced in the signal chain before the injection of phase calibration signals or due to spurious phase calibration signals (see, Rogers, 1991). As we have stated, a large fraction of the systematic differences between group delay and phase delay observables is attributed to group delay observables. It is important to identify the causes of these systematic differences, as group delays are the basic observables of geodetic VLBI. However, any conclusion on this effect needs further investigation.
(7) Polarization-related effects. The polarization leakage (Martí-Vidal et al., 2021) can have different impacts between the observations of the legacy S/X and the mixed modes. The visibility of RH + RV from the mixed mode session is not able to minimize the impact of the D term as is the pseudo-Stoke I visibility constructed for the dual linear polarizations in VGOS; without being calibrated, there may be significant errors. However, the short baselines with more than two sessions available were observed in the legacy S/X mode. The results so far do not allow us to assess the potential errors in the mixed mode of the circularly polarized receivers and the linearly polarized receivers. The only information that we presently have is that the WRMS residual level of the mixed mode observables of baselines ONSALA60-ONSA13NE and ONSALA60-ONSA13SW is about 5.0 ps and that of baseline ONSA13NE-ONSA13SW, both of which observed with the linearly polarized receivers, is only 2.0 ps based on the data analyses of the two mixed mode sessions. However, the former two baselines are about five times longer, which prevents us from drawing any conclusion. We leave these effects for a later investigation.
(8) Source positions. The short baselines in this study are sensitive to source structure only at the sub-arcsecond level or larger, three orders of magnitude greater than the typical uncertainties of source positions in the third realization of the International Celestial Reference Frame (ICRF3; Charlot et al., 2020). In general, source positions are not estimated in data analysis of short-baseline observations, as one would keep the number of model parameters to a minimum to obtain robust solutions. Due to radio interferometry, source structure at the arc-second scales is resolved out and not detected by the long baselines of thousands of kilometers. The source positions in the ICRF3, determined primarily from long-baseline observations, were used as a priori and not estimated in our solutions. However, the ICRF3 source positions can be different from the reference positions for the short baselines due to large scale structure for some of the sources. Even though we have stated that these effects should be small for most of the geodetic sources, the magnitude should be properly quantified by actual observations. We analyzed the VLBI sessions of station position tie carried out at the Onsala site, the ONTIE sessions (Varenius et al., 2021). The advantages of these ONTIE sessions for this particular investigation are the large number of scans (∼1,200 scans per session) and the multiple short baselines formed by the three antennas at the site; they allow us to obtain robust solutions for estimating a large number of source position parameters. For instance, session 20NOV12VB observed 126 sources with 3,522 observations in total. With the parameterization of one-hour-interval PWL functions for both the dZWDs and the clocks of antennas ONSA13NE and ONSA13SW and the positions of the two antennas, the WRMS delay residual is 2.9 ps based on 3,442 used phase delays. The formal errors of station position estimates are 0.2 mm in the up direction and 0.05 mm in the horizontal directions. In another solution, we estimated right ascension and declination for the 105 sources with more than 10 phase delays together with the same parameters in the previous solution. The WRMS delay residual is 2.7 ps from the solution of estimating source positions; however, the formal errors of station position estimates increase by 50% in the up direction and by 100% in the horizontal directions. The differences in station position estimates are about 0.1 mm in the up direction and 0.2 mm in the two horizontal directions, a demonstration that errors in source positions affect the horizontal directions rather than the up direction for the ties. The results from other ONTIE sessions are similar to that of session 20NOV12VB presented here. Radio sources observed in the ONTIE sessions are from the same source catalog of the IVS observations; thus, it is reasonable to conclude that the impact of source position differences due to large scale structure may be about 0.2 mm on the horizontal directions. This is insignificant relative to the uncertainties and probably will not cause systematic changes in position estimates but only increase the scatters. There is no intention either to observe the same set of sources or to observe a given source at the same Greenwich mean sidereal time in different sessions, therefore, this error source is non-repeatable (i.e., a random error as opposed to a systematic error).
(9) Atmospheric effects. In order to separate the elevation-dependent effects of the atmosphere from the estimates of station positions, in particular the Up component, geodetic VLBI observations rapidly switch between radio sources at different directions. Fortunately, atmospheric effects are greatly canceled for short baselines and the hydrostatic part of the effects is modeled in our solutions. The residual effects are investigated and discussed here.
Using phase delays of WETTZ13N-WETTZELL in the 107 sessions, we obtained three sets of solutions: (a) estimating dZWDs with a time interval of 1 hr for a PWL function, (b) estimating dZWDs with a time interval of 24 hr, and (c) not estimating dZWDs. Differences in the mean estimates of the horizontal components are negligible, smaller than 0.01 mm, between these three sets of solutions, and the mean of the differences in the Up component is 0.05 ± 0.05 mm between estimating dZWDs with two different time intervals. However, the mean of the differences in the Up estimates of antenna WETTZ13N is 1.1 ± 0.1 mm (estimating minus not estimating dZWDs). The residual time series of the Up estimates from the three sets of solutions are shown in Figure 8 and are available as a data set in Supporting Information S2. Seasonal variations with a magnitude larger than 1 mm occur in the Up residuals when the dZWDs are estimated. Even though the atmospheric turbulence due to wet troposphere on the local scale is believed to have detectable effects in VLBI observables (see, e.g., Treuhaft & Lanyi, 1987), our results suggest that the baseline time series can be stabilized by not estimating the atmospheric effects in the case of the short baseline WETTZ13N-WETTZELL at the risk of biasing the Up component. Through the investigation, we conclude that the impact of atmospheric effects on the baseline ties from short-baseline observations is negligible in the horizontal components but can be 1-2 mm in the Up component. This amount of impact may be inevitable and prevalent in the current VLBI products because the signal propagation delays due to water vapor are always estimated in geodetic solutions of VLBI observations. The study may suggest that in order to achieve global geodetic accuracy of 1 mm with VGOS alternatives for measuring wet delays may need to be further developed and verified, such as collocated microwave radiometers (see, e.g., Forkman et al., 2021).

Potential Applications of Phase Delays From Long Baselines
We have demonstrated that phase delays can be used (a) to investigate the systematic errors of group delays and (b) to derive geodetic results of baseline vectors. However, the study was limited to utilizing these observables from short baselines, mainly due to the challenge in resolving phase ambiguities for long baselines. In this section, we will briefly discuss the potential application of phase delays from long baselines.
As discussed in Section 2.2 in detail, the methods of resolving phase ambiguities include: (a) directly employing group delays to predict the phase ambiguities, (b) using theoretical model delays, and (c) relying on geodetic solutions of group delays. Because instrumental effects typically can introduce large systematic errors in the ionospheric corrections determined from S/X-band group delays (e.g., constant offsets) and resolving phase ambiguities for long baselines requires a special attention to the difference in the reference frequencies of these corrections between group delays and phase delays (e.g., the constant offsets can introduce large variations due to the changes in the reference frequencies of the corrections from group delays), none of these three options could work well for long baselines in the legacy VLBI observations. In the VGOS observations, however, because the ionospheric effects are simultaneously estimated with group delays and visibility phases in the broadband fringe fitting process (Cappallo, 2014), this naturally resolves the challenge of predicting phase ambiguities caused by the ionospheric effects. The VGOS phase delays are recently discussed and investigated by the IVS VGOS Technical Committee, and their potential applications can be to improve the quality of the group delays and to investigate the effects of source structure (Xu et al., 2022) with the goal of improving the geodetic products from VGOS. We remark that without mitigating the systematic errors caused by, for example, water vapor and source structure, the higher precision of phase delays alone cannot promise significantly better geodetic products in general cases of long baselines.

Summary
We have analyzed the phase delays in the IVS routine global observations for the short baselines at eight geodetic VLBI sites to derive the station position ties with high accuracy. The results of the baseline vector WETTZ13N-WETTZELL have baseline repeatabilities of better than 1 mm in all the three directions. The potential systematic errors were investigated and discussed in the study. As demonstrated by the investigation of the cable delay corrections, instrumental effects typically can introduce errors with a magnitude larger than 1 mm in station position estimates. The atmospheric effects, which are always estimated in geodetic solutions, may also cause seasonal fluctuations at the a few mm level in station position time series.
Phase delays produce significantly better determinations of the baseline vectors than the linearly combined group delays at S/X bands. An independent solution of only short-baseline observables does not suffer from some of the errors in long-baseline observables. The phase delay results of the position ties can be directly used in the data analysis of legacy S/X observations and VGOS observations or the combined analysis of both observations. Nevertheless, it should be noted that there currently exists incompatibility in applying these phase delay results to the routine data analysis based on group delays since some antenna pairs have significantly different position ties between group delays and phase delays. When the VLBI phase delay results are used for studies with other space techniques (see, e.g., Ning et al., 2015;Glaser et al., 2019), attention should be paid to take the systematic, repeatable errors into account, in particular the Up coordinate.

Data Availability Statement
The databases and the original visibility data of the VLBI observations used in this work are publicly available in the three primary data centers of the International VLBI Service for Geodesy and Astrometry. The repositories can be accessed at https://ivscc.gsfc.nasa.gov/productsdata/data.html. Data analysis software that was used in this work, νSolve, is available for public access at https://sourceforge.net/projects/nusolve/. The result data are provided in the Supporting Information S2.