North America presently exhibits vertical and horizontal crustal motion due to glacial isostatic adjustment (GIA). To explore these motions across central and eastern North America, the Earth model response was determined for 500 3-layered randomly generated mantle viscosity profiles, assuming the density and elastic parameter structure of the spherically symmetric Preliminary Reference Earth Model (PREM). The ICE6G_C surface loading model was utilized and motions from hydrological loading and current ice loss were accounted for. The predictions were compared to observed velocity fields downloaded from Nevada Geodetic Laboratory for selected Global Navigation Satellite System (GNSS) sites across North America. The fit is assessed through a Root Mean Square (RMS) calculation of the residual velocities. Scanning across lithospheric thicknesses, the best fitting models were compiled to assess the overall best fit for vertical, horizontal, and combined crustal motions. The horizontal and combined response exhibit two optimal viscosity profile ranges that depend on the lithosphere thickness. The profile at thinner lithosphere (<120 km) is akin to profiles like VM1 but inferred best-fitting mantle viscosities shift by an order of magnitude at thicker lithospheres (>120 km). The optimal viscosity for the vertical velocities is relatively consistent with profiles similar to VM5a and VM1. Despite the large number of tested Earth models, the residues of best-fitting models are significant and may indicate the need for more complex Earth models incorporating lateral heterogeneity to improve the fit to the observations.