Abstract: Consider an expensive computer simulation of a complex system whose inputs are governed by a known distribution and whose output demarcates passing and failing. For example, we are motivated by a simulation of airflow around an aircraft wing with inputs specifying wing design/flight conditions and output indicating whether aerodynamic efficiency standards have been met. Our ultimate objective is to quantify the probability of system failure (which could be rare) with only several hundred evaluations of the expensive simulation. We tackle this problem in three parts. First, we develop a Bayesian deep Gaussian process (DGP) surrogate to furnish predictions and uncertainty quantification at unobserved inputs. DGPs outperform ordinary GPs when dynamics are nonstationary. Second, we propose a contour locating sequential design scheme to train the DGP to identify the failure contour in the response surface. Third, we incorporate a hybrid Monte Carlo estimator of the failure probability which combines DGP surrogate predictions with strategically allocated evaluations of the expensive model. All methods are supported by publicly available software and the “deepgp” R package on CRAN.