Sometimes a treatment, such as receiving a high school diploma, is assigned to students if their scores on two inputs (e.g., math and English test scores) are above established cutoffs. This forms a multidimensional regression discontinuity design (RDD) to analyze the effect of the educational treatment where there are two running variables instead of one. Present methods for estimating such designs either collapse the two running variables into a single running variable, estimate two separate one-dimensional RDDs, or jointly model the entire response surface. The first two approaches may lose valuable information, while the third approach can be very sensitive to model misspecification. We examine an alternative approach, developed in the context of geographic RDDs, which uses Gaussian processes to flexibly model the response surfaces and estimate the impact of treatment along the full range of students that were on the margin of receiving treatment. We demonstrate theoretically, in simulation, and in an applied example, that this approach has several advantages over current approaches, including over another nonparametric surface response method. In particular, using Gaussian process regression in two-dimensional RDDs shows strong coverage and standard error estimation, and allows for easy examination of treatment effect variation for students with different patterns of running variables and outcomes. As these nonparametric approaches are new in education-specific RDDs, we also provide an R package for users to estimate treatment effects using Gaussian process regression.