Combining VGOS, legacy S/X and GNSS for the determination of UT1

We perform a combination on the observation level (COL) between VLBI and co-located GNSS in the context of VLBI intensive sessions. Our approach revolves around an estimation procedure which uses 3 h of GNSS data that uniformly encapsulate the 1-h VLBI data, in order to provide consistent troposphere information. We test this approach on both VGOS and Legacy S/X using the VGOS-B and concurrently observed INT1 sessions. The COL strategy is found to increase the precision by 15 % over both session types and leads to an increase of 65 % in the agreement between the sessions when estimating tropospheric gradients every 3 h. A more frequent estimation of the gradients every 1 h, which can be rigorously pursued with the utilization of multi-GNSS, results in a further convergence of the two session types by 30 %. The COL-aided length-of-day (LOD) products also show a 55 % better agreement to external GNSS-derived LOD. In the light of the increasing precision of broadcast GNSS orbits and clocks, this COL strategy can be used to derive rapid UT1-UTC products.


Introduction
Geodetic Very Long Baseline Interferometry (VLBI) is an interferometric technique that relies on observing a catalogue of well-defined distant radio sources (Sovers et al. 1998). These observations provide a quasi-inertial external reference which allows for the determination of the parameters that are needed for the transformation between terrestrial (earth-fixed) and celestial (inertial) reference frames. Among the full set of those parameters, which are called the Earth Orientation Parameters (EOP), the Universal Time (UT1) refers to the Earth's angle of rotation about the Celestial Intermediate Pole (CIP) axis. The rate of the rotation varies at different timescales, reflecting the effect of underlying geophysical phenomena, like ice-water mass redistribution, tectonic plate movement and planetary and lunisolar gravitational couplings, etc. (Ferrandiz et al. 2015). A precise determination of UT1 is thus essential for the separation and study of the physical processes that influence it, and it is critical in the fields that rely on accurate determination of the position and velocity of artificial satellites, like space navigation (Max-Moerbeck et al. 2015), satellite altimetry (Zelensky et al. 2018) and geodesy (Montenbruck and Steigenberger 2020).
For this reason, and with the aim of studying highfrequency UT1 fluctuations, the International VLBI Service for Geodesy and Astrometry (IVS) (Nothnagel et al. 2017 regularly organizes the so-called intensive (INT1,INT2,INT3) sessions. These are 1-hr-long observing sessions that routinely involve two (INT1) or three (INT2 or INT3) stations which form long east-west oriented baselines. This particular geometry enables high-accuracy determination of short-term variations of UT1, which are reflected on the estimation of (UT1-UTC), the difference between the UT1 and the coordinated universal time (UTC). While long east-west baselines have been accepted as providing the highest sensi-tivity in (UT1-UTC) estimation, recent developments show that a different inclination of the baselines with respect to the equatorial plane may be needed to achieve enhanced accuracy (Schartner et al. 2021).
The development and integration of the next-generation VLBI Global Observing System (VGOS) antennas into the IVS observation network (Behrend et al. 2019), motivated new types of observing sessions to be organized to take advantage of the unique features of VGOS (Petrachenko et al. 2009). The broadband receivers and the relatively smaller antenna size with respect to the legacy S/X telescopes, enable shorter integration times and faster slewing manoeuvres which translate to a bigger number of observations and a more comprehensive and geometrically diverse sky coverage. In this context, intensive-like sessions were performed between December 2019 and March 2020 that did not involve the legacy S/X VLBI network but instead utilized VGOS stations at Onsala (Sweden) and Ishioka (Japan). The so-called VGOS-B sessions were scheduled to be concurrent to regular INT1 sessions, using a scheduling approach that favoured the observation of sources at the corners of the mutually visible sky among the stations . The coordination between VGOS-B and INT1 sessions, and their purposeful scheduling design which capitalizes on the features of the respective observing systems to provide an accurate (UT1-UTC) estimation, makes a comparison between them meaningful as it illustrates the expected improvements that one can anticipate from VGOS. In fact, such a comparison has already been made and showed that the VGOS-B sessions were 30-40 % more accurate in terms of root-mean-squared error and standard deviation , reflecting the results of previous simulations on the topic (Corbin and Haas 2019).
We intend to broaden the comparison scope between the VGOS-B and concurrent INT1 sessions, introducing various inter-and intra-technique combination schemes. In particular, we utilize the combination on the observation level (COL) approach to estimate common tropospheric parameters for the co-located VLBI and Global Navigation Satellite Systems (GNSS) stations, using the c5++ multi-technique space-geodetic analysis software .
The limited number of observations of the intensive sessions only allows for the zenith wet delay to be modelled as a constant offset while no estimation of tropospheric gradients is possible. We thus expect that the improved resolution of the tropospheric parameters that the COL will offer, should positively affect the determination of (UT1-UTC), either through direct correlation or indirect one, e.g. by a more precise estimation of the clock bias. We also combine the concurrent INT1 and VGOS-B sessions on the basis of a common (UT1-UTC) estimation, essentially integrating all baselines in one solution, with the aim of further increasing the sensitivity.
Section 2 provides a description of the data used and presents the estimation setup with the parameterization of the different combination schemes. This is followed by Sect. 3, which presents the acquired accuracy for the different combination options, Sect. 4 which provides a discussion of the results with the acquired insights, and, finally, Sect. 5 with a summary of the key findings and an outlook for further research.

Methods
We analyzed 11 out of the 12 total VGOS-B sessions and their concurrent INT1 sessions between December 10, 2019, and February 24, 2020. The session on January 7th was significantly degraded due to the fact that VGOS bands B/C/D were missing in ISHIOKA, and thus was excluded from the analysis.

Input data
The data were acquired from the CDDIS data centre (Noll 2010) in the so-called vgosDB format (Gipson 2014). This is a netCDF compliant database which splits sessiondependent, station-dependent, station/scan-dependent, and observation-dependent data in separate bins for efficacious modular processing. We used databases of version-4 or higher which contain resolved group-delay ambiguities, and in the case of Legacy S/X, ionosphere corrections. The latter is not required for data from the VGOS antennas, since their increased bandwidth allows for the compensation of ionospheric-related effects in the fringe-fitting process, before the generation of the vgosDB databases. The codes of the 11 analyzed VGOS-B and INT1 observing sessions, along with the participating stations, are shown in Table 1.
We also analyzed data from the co-located GNSS stations, using the precise point positioning (PPP) technique. The data were attained in the so-called Rinex-3 format (Gurtner and Estey 2007), a receiver-independent format which enables multi-GNSS processing (Montenbruck et al. 2017). For the purposes of this study, we utilized our internally developed multi-GNSS compliant module that can be used in tandem with the c5++ analysis software and expands its GPS-only processing capability. The timing inconsistencies among the GNSS constellations, known as inter-system biases, are internally resolved within the estimation process to ensure homogenization of the multi-GNSS observables. The clock and orbit products that are necessary for the PPP processing were obtained from the International GNSS Service (IGS) (Johnston et al 2017). In particular, we used the precise orbits and clocks produced by the Center for Orbit Determination in Europe (CODE) (Prange et al. 2015). The participating GNSS stations and constellations processed are

Analysis options and parameterization
The analysis options can be separated into three categories, namely (a) single-network solutions, (b) single-network COL with GNSS and (c) dual-network COL with GNSS. The general solution codes along with the parameterization are presented in Table 3. We performed several variations of the solutions which will be further denoted with an underscore on the general solution code. Two aspects were taken into account when designing the approach to the GNSS-aided COL. The first was the starting time and the duration of the analysis. During conventional processing of a VLBI intensive session, a priori (UT1-UTC) is obtained via Lagrangian interpolation at the middle time node. For consistency reasons we sought to make sure that the middle time node of the COL coincides with that of the 1-h INT1 and VGOS-B sessions. This means that, for every session, the analysis window should expand equidistantly from 19 UT on each side. We wanted to keep the "intensive" concept, i.e. a short time interval with high sensitivity to the estimated parameters, and expanded with 1.5 hr on each side. We essentially created a 3-hr GNSS "shell" which encapsulates uniformly the 1-hr VLBI session. In this way, the tropospheric parameters of interest that apply to this 1h VLBI session are not boundary parameters, ensuring the consistency of the estimates.
The short time window also allowed us to study the impact of multi-GNSS in greater detail. In an environment of limited contributions from each GNSS constellation observationwise, we expect that the role of geometry will become the dominant factor which will influence the results. This multi-GNSS aspect of the analysis was the second major factor when designing the framework, and its effect is presented in detail in Sect. 3.
The station coordinates were kept fixed in the analysis, except for the stations ISHI and ONS1 where no ITRF2014 coordinates were available. The d(UT1-UTC) estimate was treated as a constant offset, and refers to the correction that was estimated with respect to the a priori EOP. The clock offset was estimated as a 2nd-degree polynomial, in the case of VLBI, and piecewise-linear (PWL) offsets, in the case of GNSS. For one session, the MKEA clock derived from a 24-h PPP solution was given as a priori, since it was found to exhibit significant short-term variability that could not be sufficiently resolved with 3 h of GNSS data. The parameterization was the same as for the 3-h sessions with the exception of solving for station coordinates as well. The zenith wet delay (ZWD) was estimated as PWL offsets, while the east (EG) and north (NG) gradients were treated as constant offsets (CO) when combining VLBI with a single GNSS constellation, but estimated as PWL offsets for combinations with two or more GNSS constellations. The zenith hydrostatic delay (ZHD) was modelled using the Saastamoinen formula (Saastamoinen 1972), and the ZWD bias (δZWD) between co-located stations was separately computed and included as a priori information according to established practice (Rothacher et al. 2011). Rate constraints were applied in the clock offsets and ZWD. The details on the general parameterization setup are presented in Table 3.
We used a batch least-squares (LS) estimator with variance component estimation (VCE) to ensure the elimination of discrepancies between the employed noise model and the error residuals and consequently to weight in a rigorous manner the relative contribution of the different observables. The observation bins that the VCE was based on were, (a) per observation type, (b) per GNSS station, (c) per GNSS constellation and (d) per phase or code GNSS observable. The GNSS observations were downsampled to the 5 min resolution. The IVS VTRF2020b (BKG 2020) and the ITRF2014 (Altamimi et al. 2016) were used as a priori values for the VLBI and GNSS station coordinates, respectively. The VTRF2020b coordinates that were determined through dedicated shortbaseline interferometry campaigns, were used for the VGOS Onsala twin telescopes ). The ICRF3 (Charlot et al. 2020) realization of the celestial reference frame and the IERS-14-CO4 series (Bizouard et al. 2019) was used as a priori information on source positions and EOP. Displacements that occur due to ocean, solid and pole tides (Petit and Luzum 2010), S1-S2 atmospheric pressure loading (Ray and Ponte 2003) and post-seismic deformations (Altamimi et al. 2016), were included in the modelling.

Tropospheric parameters during COL and intensive sessions
In order for the COL strategy to be meaningful, the VLBI and GNSS co-located stations must estimate consistent tropospheric parameters. Those are derived from the refractivity of the troposphere which is split into the (stochastic) "wet" and (well-modelled) "dry" parts, the latter being dependent on latitude, ellipsoidal height and local pressure. Consistent inter-technique meteorological data are thus essential for the COL to be successful. While local weather data are available for VLBI as they are included in the vgosDB database, the c5++ analysis software provides consistency through the use of the Global Pressure and Temperature (GPT3) model (Landskron and Böhm 2018) across all space-geodetic techniques.

On the possibility of intra-VGOS troposphere ties
We estimated common ZWD, NG and EG between colocated GNSS (G) and VLBI (V) stations modelling the slant delays (SD) as where m f w , m f d and m f g are the wet, dry and gradient mapping functions, and α is the azimuth angle. The VGOS-B sessions offered an option for such a combination between the two different VGOS antennas, Oe and Ow, at the Onsala site. However, the common tropospheric parameters coupled with the observations of the baseline that those two antennas formed, would introduce an ambiguity in the troposphere definition. To illustrate this, we consider the scenario where a single ZWD for both Ow and Oe is estimated. After removing the gradients, and the modelled ZHD and δZWD, the residual delay due to troposphere for the Oe-Ow baseline would be where δm f w,Oe−Ow is the now-modified mapping function and a source of inconsistency. It is unbounded, contrary to mapping functions that exhibit a lower bound at 1 as to reflect the (minimum) refractivity at zenith. It is dependent not only on the elevation angle, , but also on the difference of the elevations of the co-located stations, δ . It takes values very close to zero, when the two stations comprising the baseline are in close proximity, which makes it a poor predictor of incremental ZWD changes. To avoid such ambiguities in the formulation of the estimation problem, one could either remove the Oe-Ow baseline from the data processing and combine Oe, Ow with co-located GNSS or estimate separate tropospheric parameters for the VGOS stations and combine them with different co-located GNSS antennas. We chose to do the latter, as we found out that the single-technique solution comprising of all available baselines was the most precise one in terms of achieved formal errors, and subsequently chose to use that as our reference solution.

Tropospheric gradients
The tropospheric gradients, and in particular the one that lies on the east-west direction, have been shown to have a significant impact in the quality of (UT1-UTC) (Böhm and Schuh 2007) and can amount to some tens of μs (Böhm et al. 2010). In fact, previous research points towards including a priori determined gradients, and estimating ZWD, as the preferred analysis scheme for VLBI intensive sessions (Teke et al. 2015), and the effect of ray-traced gradients on (UT1-UTC) has also been explored (Landskron and Böhm 2019).
In this study, we relied on the COL to deliver consistent tropospheric parameters, and expect that the explicit presence of the correlations in the variance-covariance matrix would further improve the accuracy of the (UT1-UTC) estimates.
We investigated the effect of tropospheric gradients by altering their parameterization between 3-hr and 1-hr PWL. We also compared our approach of analysing VLBI intensive sessions by means of an encapsulating COL with GNSS, to the state-of-the-art approach of including a priori computed gradients.

Additional data processing for length-of-day estimation
We derived length-of-day (LOD) estimates using that enabled a comparison to external LOD products. This calculation necessitated the COL to be performed for a secondary dataset that of the day which follows each day of the concurrent VGOS-B and INT1 sessions. As no VGOS-B sessions existed then, the LOD estimation was performed using the INT1 sessions. In two cases, I19357 and I19364, there were no follow-up sessions the next day, and consequently they were excluded from the LOD computation. The secondary dataset is presented in Table 4.

Results
The time series of the initial single-technique d(UT1-UTC) estimation for INT1 and VGOS-B is presented in Fig. 1. Three main statistics are visible, namely the weighted mean offset (WMO), which is a weighted mean taken from the

Inter-technique combination
We first examine the inter-technique combination results for the INT1 and VGOS-B sessions using 3-hr tropospheric gradients. The effect of incrementally adding the different GNSS constellations to the V1 solution, in terms of STD and MFE, is presented in Fig. 2. Notice that, in this case, the generic solution code VG1 is supplemented by a subscript to denote the number of constellations processed. In particular VG1 1c means that only GPS is utilized, VG1 2c refers to GPS and Galileo, VG1 3c contains GPS, Galileo and GLONASS, while VG1 3c+ all available data. The generic solution code, VG1, implicitly refers to VG1 3c+ , and the same happens in an equivalent fashion for VG2. The effect of stepwise addition of GNSS constellations to solution V2 becomes clearly visible when one estimates tropospheric gradients more frequently than the default 3-h resolution. Hence the relevant plot is presented in Sect. 3.3. The effect of COL is studied with respect to derived (UT1-UTC) both in precision and accuracy. While the former can be evaluated using the STD and MFE, the latter is harder to assess, given that the INT1 sessions have already been incorporated in the generation of the IERS final products (Luzum and Gambis 2014).
Our approach to gauge the accuracy of our solution is twofold. Firstly, we rely on the fact that the INT1 and VGOS-B sessions are two independent observation sets, and their COL is performed using the respective co-located GNSS stations. This independency can be the basis of the generation and evaluation of an accuracy metric. Before any intra-technique combination is attempted, which functionally correlates these two sessions on the basis of a common normal equation, one should compare whether the intertechnique combination results in the d(UT1-UTC) estimates  Fig. 3 in the weighted root-mean-square (WRMS) differences. The same effect can be seen but not readily quantified in Fig. 1, where the degree of similarity of the obtained WMO can also be examined. The effect of COL on the accuracy can also be determined using external products. The GNSS-derived length-of-day (LOD) estimates from two analysis centres, IGS (Kouba 2015) and CODE (Dach et al. 2016), are used as the basis of the comparison. In order to obtain LOD estimates from our analysis setup, we additionally perform the COL and d(UT1-UTC) estimation for the day (t 2 = t 1 +1 d) that follows every day (t 1 ) in our dataset. We then compute LOD at the middle time node t = t 2 −t 1 2 using Eq. 1. The INT1 sessions typically provide a (UT1-UTC) estimate at 19 UT, which means that the processing of subsequent days in pairs results in an LOD estimate at 7 UT. This comes in contrast to the GNSSderived daily LOD products that are given at 12 UT. We use a 5-point Lagrangian interpolation to resolve this discrepancy and temporally align the two LODs. The standard deviation and the mean bias are then calculated between the solution sets, and the results are presented in Table 5.
The incremental effect of the addition of the GNSS constellations in the estimation of the gradients is given in Table 6. More specifically the mean formal error and mean bias of the estimates are given for solutions VG1 and VG2.   The impact of each separate gradient component on d(UT1-UTC) is presented in Figs. 4, 5, where the differences in the d(UT1-UTC) estimates between single-technique solutions (V1,V2) and combined solutions with all constellations (VG1,VG2) are plotted against the sum of the total east and north gradients over the stations. We also include the 9 additional sessions that were used to extract the LOD estimates.

Intra-technique combination
The combination of VG1 and VG2 solutions on the basis of a common d(UT1-UTC) estimate provides an augmented solution denoted as VG12. The intra-technique combination results are presented in Fig. 6a, where the time series of the original VG1 and VG2 solutions is plotted together with that of the VG12. The WRMS differences of VG12 with respect to the two inter-technique COLs are presented in Fig. 6b. (b) The WRMS differences (inµs) of the solution VG12, with respect to VG1 and VG2, as denoted in the y axis. Fig. 6 The effect of the intra-technique COL of VG1 and VG2 (VG12) in the time series and WRMS differences

Estimation of 1-h gradients
We take advantage of the improved sky coverage that multi-GNSS gives to study the asymmetricity of the troposphere in more detail. In particular, we are interested to see if denser resolution of the tropospheric gradients affects the quality of our estimation. We generate two modified solutions for VG1 and VG2, denoted as VG1 1h and VG2 1h , with increased resolution for the PWL model of the gradients from 3 hrs to 1 hr. The results in the form of time series and WRMS differences are presented in Fig. 7. The increased precision of the estimated parameters, enables us to investigate the effect that the inclusion of more GNSS constellations has in the VG2 solution, and is shown in Fig. 8. Notice that the STD metric of Fig. 2

Comparison to single-technique analyses with a priori gradients
A single-technique solution where tropospheric gradients are provided as a priori information and ZWD is estimated online, has been described as the most precise formulation for the processing of intensive sessions (Teke et al. 2015). We validate this assumption and evaluate the performance of COL against it. We first analyze the 24-h GNSS data for the corresponding days and stations, and estimate ZWD and gradients. Two modified V1 and V2 solutions are performed where the first contains the a priori computed gradients, and the second, in addition, the a priori estimated ZWD. We then use VG1 and VG2 as the bases for a cross-comparison. More specifically, to evaluate the modified V1 together with the VG1, we compare their WRMS differences with respect to VG2. In a mirroring sense, to evaluate the modified V2 together with the VG2, we compare their WRMS differences with respect to VG1. The results of this comparison are presented in Fig. 9. Notice that the suffix g stands for gradients and t for the whole troposphere parameter set. (b) The WRMS differences of VG2, V2, and its augmentations with gradients (V2g ) and the complete troposphere set (V2t) with respect to the VG13c+ solution. Fig. 9 Comparison between COL and single-technique solutions that are augmented with externally estimated troposphere information

Discussion
The VGOS-B sessions achieved a twofold increase in observations per baseline compared to the INT1 sessions and with significant attention to using sources at the edges of the mutually visible sky . The results of this purposeful scheduling design can be seen in Fig. 1a. VGOS shows a 78 % decrease in the MFE of the estimates, compared to the legacy S/X. This is a marked increase in the sensitivity of the measurements to incremental changes in (UT1-UTC) and positively affects the accuracy of the estimates. We rigorously formulate the argument as to why this is the case later in the discussion. The application of the COL, as seen in Fig. 1b, causes a decrease of approximately 15 % in the MFE across both sessions. Other effects are the convergence of the WMO with respect to each other and also towards the zero level which indicates that the COL is inducing an alignment effect to UT. The scatter of the estimates, as tracked by the STD metric, is also significantly reduced. The INT1 sessions are shown to be more affected with a 48 % decrease, while the VGOS-B sessions show a 15 % decrease in STD. What can be extracted from Fig. 1 is that the COL provides a more precise estimation of (UT1-UTC). We want to see whether this statement can be extrapolated to include the accuracy of the estimates. We first inspect an increased correlation of the two time series in the COL solution. The degree of this correlation is given in Fig. 3, in the form of WRMS differences. In particular, the COL causes a 62 % decrease in the WRMS between the two types of sessions. The biggest modulating effect is exerted on the (UT1-UTC) estimates of the INT1 sessions where the combined solution shows a difference of 21.3 μs with respect to the single-technique solution. The COL on the VGOS-B sessions results in a difference of 6 μs with respect to the single-technique solution. The COL is therefore causing these two independent sessions to converge with respect to each other. It can thus be extracted that the combination is increasing the accuracy of the estimation. This is further consolidated when an evaluation with an external reference is attempted, as it is given in Table 5. The comparison to CODE and IGS products shows that the COL is causing a 58 % and 51 % increase in the agreement, respectively. The result is that the STD differences fall in the 20 μs level. The mean bias of the estimated LOD is on the sub-μs level with respect to the IGS product, while it is around 5 μs for CODE. This may be explained by the fact that the CODE product shows a systematic bias as it is computed solely from GNSS observations (Zajdel et al. 2020). On the other hand, the STD metric shows a better agreement to CODE which might reflect its enhanced precision. Overall, the COL is increasing the accuracy of the estimates of both INT1 and VGOS-B. A smaller normalized effect is seen in VGOS-B which can retrospectively indicate that it is inherently more accurate.
The incremental addition of the GNSS constellations stabilizes the gradient estimates and lowers their formal errors. This can be clearly seen in Table 6, where for both VGOS-B and INT1 the mean bias converges with the incorporation of 2 or more constellations. The formal uncertainties get 40 % lower in the case when all constellations are used as opposed to GPS-only, and fall to the 0.2 mm level. Thus, using a multi-GNSS estimation scheme allows for incremental changes in gradients to be more visible and robustly determined. The importance of this is reiterated in Fig. 4, where the high correlation of d(UT1-UTC) to the east gradient is presented. Consistent east gradient estimates are essential to a precise d(UT1-UTC) determination. The north gradient exhibits little correlation to d(UT1-UTC) as shown in Fig. 5 which coincides with findings of previous studies (Teke et al. 2015). The effect that multi-GNSS processing has in the estimation can be seen in Figs. 2 and 8. The results show that adding more constellations progressively lowers the MFE of the estimates. The INT1 sessions are more sensitive to this inclusion, as an improvement can be already seen when estimating 3-hr gradients. This is consistent with the fact that the correlation of d(UT1-UTC) and tropospheric gradients is well established (Böhm and Schuh 2007;Teke et al. 2015), and it is stronger in the INT1 sessions due to the lower amount of overall observations and hence the decreased over-determination of the estimation problem. The improvement of the MFE and STD can be seen in Fig. 2 and follows an asymptotic pattern with more constellations being added. This is both expected as tropospheric gradients can be determined with a finite precision, but also an implication of the fact that not all antennas detect BDS and QZSS, as can be seen in Table 2. This in turn leads to them having a decreased share in the total amount of observations. The effect of incrementally adding constellations in VGOS-B was a lot less discernible for the analyses with 3-h gradients. The biggest gain happened when the gradients were first introduced in the estimation, but further refinement showed diminishing returns. In order to get the asymptotic nature seen in Fig. 8, we switched to a denser 1-hr gradient estimation. It seems to follow that the VGOS-B sessions, as intrinsically more accurate that the INT1 ones, need more accurate gradient estimation in order for significant improvement to be detected. In this context, we also swapped the less precise STD metric for the WRMS differences with respect to the, now well established as reasonably accurate, VG1 solution. Notice that the WRMS in this case falls from the initial 12.2 μs for V2 to 8.2 μs for VG2 3c+ , which is slightly increased with respect to the lowest detected value of 8.1 μs that was obtained with VG2 3c . In general, more precise gradients showed a positive effect on both sessions. This can be clearly seen both in terms of time series in Fig. 7a and in terms of WRMS where the difference between the VG1 and VG2 solutions decreases by 30 % when switching from 3-hr to 1-hr gradients, as seen in Fig 7b. The effect of esti-mating gradients with higher temporal resolution has been previously studied in terms of, e.g. station position repeatabilities (Zhou et al. 2017;Lu et al. 2016) where improvements of up to 17 % have been found. We consider our results to be in line with these improvements, given the limited dataset used in the present study. We believe that the results shown motivate further studies in high-frequency tropospheric gradient modelling and estimation.
We already saw the effect of the inter-technique COL in alignment of the estimates to UT. An intra-technique combination of both legacy S/X and VGOS networks holds the potential to determine the d(UT1-UTC) time series in the most consistent way, as in a sense this combination scheme provides a weighted mean of the results of the two constituent networks. This is indeed shown in Fig. 6a, where both by inspection of the time series and the WRMS metric, the intra-technique COL is located between the two singlenetwork solutions. It is, however, dominated by the VGOS-B estimates, as the WRMS difference to them is on the sub-1 μs level.
We have established that the COL with GNSS is beneficial for the determination of (UT1-UTC). A comparison to the current state-of-the-art approach, where gradients are given in the single-technique solution as a priori information, is presented in Fig. 9. We reaffirm previous findings that tropospheric gradients are correlated to d(UT1-UTC) estimates, and that further introducing a priori ZWD does not significantly benefit the solution. In fact, we detect a slight worsening, as Fig. 9a and 9b show. In both INT1 and VGOS-B solutions, the COL was at the same level or even outperformed the conventional approach. More specifically, the VG1 solution shows a 15 % decrease in WRMS differences to VG2, when compared with the augmented single-technique solution V1 g . In addition, the VG2 solution equivalently shows a 38 % decrease when compared to the augmented single-technique solution V2 g . The COL, as opposed to the conventional approach, retains the correlations between the estimated parameters, and this likely gives the advantage in a consistent and precise d(UT1-UTC) estimation.
The introduction of GNSS in a combination scheme with VLBI intensive sessions, like the one implemented above, aims to provide precise troposphere information in the estimation problem. Especially in terms of gradients, enhanced precision on the tenth-of-a-millimetre level is required. This means that matching the 1-hr VLBI session with GNSS data of equivalent time length may not be enough to provide the necessary consistency. One should keep in mind that, in the case of GNSS, the resolution of carrier phase ambiguities is also necessary in the context of PPP for geodetic-grade estimation. An encapsulating "shell" of 3 hrs of GNSS data and multi-GNSS processing was utilized to achieve the precision needed. This ensures the consistency of the tropospheric parameters, especially those that reside at the middle point of the (expanded) estimation problem, i.e. the ones that are present in the VLBI intensive session. Increasing the size of the "shell", need not lead to more precise results. While it may aid the determination of the tropospheric parameters themselves, it inadvertently introduces correlations between the d(UT1-UTC) parameter and all the parameters that exist outside of the 1-hr VLBI session. The choice of the size of the "shell", at least for batch estimators, is thus a matter of balance between ensuring precise troposphere in the time window of interest and limiting the number of not physically meaningful correlations that get introduced in the combination.

Conclusion and future work
In this study, we formulated a strategy for integrating the COL in the d(UT1-UTC) estimation and the intensive sessions. We utilized the independent but concurrent VGOS-B and INT1 sessions as a proving ground of the strategy. We tested different parameterization schemes and compared our results to externally derived d(UT1-UTC) products. The COL has been shown to provide precise, accurate and consistent products. The incorporation of multi-GNSS is essential in providing highly reliable tropospheric parameters, which are necessary for the COL to be meaningful. At its utmost precision, the COL matches or even outperforms the conventional augmented approach to the processing of the intensive sessions.
Although the dataset used in this study is relatively small, the improvement is consistent and tangible. To quantify it with a greater degree of certainty, an increased number of sessions should be considered. It should also be noted that the potential of the BDS constellation cannot be fully explored yet, as the ephemerides that CODE offers are limited to BDS-2 satellites. A future introduction of BDS-3 orbit and clock products will further increase the access that multi-GNSS gives to an expanded sky coverage, which will translate to a better mapping of local asymmetricities in the troposphere. This study was conducted using the precise ephemerides for the PPP solution, which inherently limits the availability of a potential (UT1-UTC) product to that of the orbit and clock products themselves. The introduction of multi-GNSS, however, holds the promise of enabling high-precision real-time PPP to be performed with the use of broadcast or ultra-rapid ephemerides (Li et al. 2015;Kazmierski et al. 2018).
We intend to further refine our data analysis incorporating LOD estimation, integer ambiguity resolution and the explicit estimation and removal of the GLONASS interfrequency code biases. We also want to try different estimation techniques, like, e.g. a Kalman filter which could give a more precise gradient estimation, or incorporating mea-surements from water vapour radiometers. We are interested in using this framework to evaluate different IGS and IVS products for a priori station coordinates, or different EOP products against the CO4 series, like e.g. the so-called finals as produced by the IERS rapid service/prediction centre (Stamatakos et al. 2020).
Author Contributions PKD designed the study, implemented the necessary extensions in the analysis software, analysed the data, interpreted the results, prepared all graphical material and wrote the manuscript. RH, EV, SM and MS organised VGOS-B series, in coordination with the IVS Coordinating Center. MS scheduled the VGOS-B experiments. EV performed the VLBI data correlation and post-correlation analysis and created the final vgosDB.
Funding Open access funding provided by Chalmers University of Technology. This work is supported by Rymdstyrelsen, the Swedish National Space Agency (Grant No. 166/16).

Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflict of interest
The authors declare that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.