A robust method to estimate vertical crustal motions by combining geocentric sea level measurements from decadal (1992 - 2003) TOPEX/POSEIDON satellite altimetry and long-term (> 40 years) relative sea level records from tide gauges using a novel Gauss-Markov stochastic adjustment model is presented. These results represent an improvement over a prior study (Kuo et al. 2004) in Fennoscandia, where the observed vertical motions are primarily attributed to the incomplete Glacial Isostatic Adjustment (GIA) in the region since the Last Glacial Maximum (LGM). The stochastic adjustment algorithm and results include a fully-populated a priori covariance matrix. The algorithm was extended to estimate vertical motion at tide gauge locations near open seas and around semi-enclosed seas and lakes. Estimation of nonlinear vertical motions, which could result from co- and postseismic deformations, has also been incorporated. The estimated uncertainties for the vertical motion solutions in coastal regions of the Baltic Sea and around the Great Lakes are in general < 0.5 mm yr^-1, which is a significant improvement over existing studies. In the Baltic Sea, the comparisons of the vertical motion solution with 10 collocated GPS radial rates and with the BIFROST GIA model show differences of 0.2 ± 0.9 and 1.6 ± 1.8 mm yr^-1, respectively. For the Great Lakes region, the comparisons with the ICE-3G model and with the relative vertical motion estimated using tide gauges only (Mainville and Craymer 2005) show differences of -0.2 ± 0.6 and -0.1 ± 0.5 mm yr^-1, respectively. The Alaskan vertical motion solutions (linear and nonlinear models) have an estimated uncertainty of ~1.2 - 1.6 mm yr^-1, which agree qualitatively with GPS velocity and tide gauge-only solutions (Larsen et al. 2003). This innovative technique could potentially provide improved estimates of the vertical motion globally where long-term tide gauge records exist.