A method of automatically combining geodetic network solutions to produce station coordinate time series with realistic computed errors has been developed and tested and is being applied on a weekly basis to Global Positioning System (GPS) global and regional networks of the International GPS Service. Our techniques include modified Helmert blocking, stochastic modeling to minimize frame bias, Monte Carlo simulation, variance component estimation, and multiparameter data snooping. An 18-month time series evaluation of 150 globally distributed stations demonstrates that our combined weekly solution is more complete, precise, and reliable than any contributing solution. Our method of attaching regional networks without perturbing the global network solution, rather than combining normal equations, improves the quality measures. The median RMS of station position residuals with respect to a constant velocity model is 2.4 mm in latitude, 3.0 mm in longitude, and 7.2 mm in height. Our solution has since been incorporated into the reference frame ITRF96 (International Terrestrial Reference Frame 1996), showing a RMS coordinate difference of 5.4 mm, the lowest of all contributing solutions. As an independent test, the RMS difference with the ITRF94 is 4.5 mm in horizontal and 8.1 mm in height. As a second external test, the station velocity solution was used to estimate plate tectonic Euler vectors, which were then compared with the NUVEL-1A model and found to differ at a level consistent with the computed errors. Given a few more years of data, our error model predicts solutions that will be sufficiently precise to rigorously test NUVEL-1A or its successors. ¿ 2000 American Geophysical Union