The harsh operating environment in a wastewater treatment process (WWTP) makes sensor faults commonplace. Detecting these faults can be challenging due to the complex process dynamics, unknown inputs, and general noise in the process and measurements. Comparing sensor readings against predictions from a physics-based or data-driven model of the WWTP is a common strategy for detecting such faults. In this work sensor measurements are directly modelled using Gaussian process (GP) regression, a data-driven multivariate approach. These GP sensor models are, with a generalised product of experts, combined into a dedicated fault isolation scheme resembling traditional observer bank methods. The residuals are monitored with a multivariate exponentially weighted moving average chart which is used for fault detection and isolation. The method is evaluated using simulated data generated with the Benchmark Simulation Model No. 1 WWTP. Fault detection performance is reported using several standard metrics such as false alarms, missed detections, time to detection, and successful fault isolations, with emphasis on reporting across a wide range of sensors and faults to provide a point of comparison for future studies. The proposed approach performs well across these metrics. Given sufficient data representative of normal operation, this approach can easily be adapted across a wide variety of plant configurations and can be used to create operatorfriendly diagnostics resembling classical control charts.