© 2019, © 2019 American Statistical Association, Institute of Mathematical Statistics, and Interface Foundation of North America. When using univariate models, goodness of fit can be assessed through many different methods, including graphical tools such as half-normal plots with a simulation envelope. This is straightforward due to the notion of ordering of a univariate sample, which can readily reveal possible outliers. In the bivariate case, however, it is often difficult to detect extreme points and verify whether a sample of residuals is a reasonable realization from a fitted model. We propose a new framework, implemented as the bivrp R package, available on CRAN. Our framework uses the same principles of the simulation envelope in a half-normal plot, but as a simulation polygon for each point in a bivariate sample. By using algorithms of convex hull construction and polygon area reduction, we describe how our method works and illustrate its functionality with examples using simulated bivariate normal data and real bivariate count data. We show how different model diagnostics can produce different results and pinpoint potential drawbacks of our approach, such as the limitations in terms of computational burden. Supplementary materials for this article are available online.