Abstract: Gaussian process (GP) surrogates are widely used for emulating expensive computer simulators, and have led to important advances in science and engineering. One challenge with fitting such surrogates is the costly generation of training data, which can require thousands of CPU hours per run. Recent promising work has investigated the integration of known boundary information within GP surrogates, which can greatly reduce its training sample size and thus its computational cost. There is, however, little work exploring the critical question of how such simulation experiments should be designed given boundary information. We propose here a new class of space-filling designs, called boundary maximin designs, for effective GP surrogate modeling with boundary information. Our designs rely on a new space-filling criterion that is derived from the asymptotic D-optimal designs of the boundary GPs from Vernon et al. (2019) and Ding et al. (2019), which can incorporate a broad class of known boundary information. To account for effect sparsity, we further propose a new boundary maximum projection design that jointly integrates boundary information and ensures good projective properties. Numerical experiments show the improved surrogate performance of boundary-integrated GPs using the proposed boundary maximin designs compared to the state-of-the-art.