Gaussian processes (GPs) are widely used as surrogate models for emulating computer code, which simulate complex physical phenomena. In many problems, additional boundary information (i.e., the behavior of the phenomena along input boundaries) is known beforehand, either from governing physics or scientific knowledge. While there has been recent work on incorporating boundary information within GPs, such models do not provide theoretical insights on improved convergence rates. To this end, we propose a new GP model, called BdryGP, for incorporating boundary information. We show that BdryGP not only has improved convergence rates over existing GP models (which do not incorporate boundaries), but is also more resistant to the "curse-of-dimensionality" in nonparametric regression. Our proofs make use of a novel connection between GP interpolation and finite-element modeling.