Elevations are one of the positional attributes embedded in all geospatial data. They are essential for a wide range of engineering and scientific activities. Some of the activities requiring precise elevations are activities of high societal impacts, such as sea level rise, storm surges and coastal inundation, floods and evacuation route planning, and crustal motion, subsidence, and other surface deformations due to seismic, mining, or other events. In order to successfully monitor and manage such events regionally (e.g., floods) or globally (e.g., sea level rise), the elevation information needs to not only be accurate but to also refer to the same zero-height reference surface (vertical datum). Accurate elevations can be obtained using spirit leveling, or by combining Global Navigation Satellite Systems (GNSS) positioning methods with a model of the geoid (the equipotential surface of the Earth's gravity field approximating the idealized mean sea surface). So although the high accuracy requirements can in general be met, the need for a common reference surface often cannot. This is because there are currently over 100 different vertical datums around the world. Their unification is therefore a scientific problem of high practical significance.Height system unification is the process of determining the offsets between existing vertical datums, each of which is defined with a fundamental zero-elevation point and the equipotential (level) surface of the Earth's gravity field that passes through that point. These offsets can be estimated if the geodetic (from GNSS), physical (from spirit leveling), and geoid (from gravity information) heights are known for the fundamental and/or other points in each datum. GNSSsurveyed leveling benchmarks on land, GNSS-surveyed tide gauge stations at the coasts, and a mean dynamic topography model at sea can be used as data, together with a geoid model. The need to estimate geoid undulations implies the solution of a geodetic boundary value problem (GBVP); see, e.g., Colombo (1980), Rummel and Teunissen (1988), Rapp and Balasubramania (1992), Venuti (2002), andMoritz (2013). Although normal heights and height anomalies from the solution of Molodensky's GBVP can and have been used for datum unification, here we work with orthometric heights and geoid undulations from the solution of the classical (Stokes) GBVP. We first review briefly the key concepts of leveling and vertical datum definition and realization, and then discuss, in spherical approximation, the single (global) vertical datum GBVP and the multiple vertical datum GBVP following primarily the developments of Rummel and Teunissen (1988). We also introduce the functional and stochastic models used for the least-squares estimation of datum offsets and their accuracies. Using an example from North America with tide gauges, we address a few important practical aspects, such as the magnitude of the so-called indirect bias term resulting from the use of biased local gravity data, and the commission and omission geoid un...