We describe the application of the phase averaging technique to the nonlinear rotating shallow water equations on the sphere, discretised using compatible finite element methods. Phase averaging consists of averaging the nonlinearity over phase shifts in the exponential of the linear wave operator. Phase averaging is a form of heterogeneous multiscale method that aims to capture the slow dynamics in a solutionthat is smoother in time (in transformed variables) so that larger timesteps may be taken. Following Peddle et al (2019), we consider finite width phase averaging windows, since the equations have a finite timescale separation. In a numerical implementation, the averaging integral is replaced by a Riemann sum, where each term can be evaluated in parallel. This creates an opportunity for parallelism in the timestepping method, which we use here to compute our solutions. Here, we focus on the stability and accuracy of thenumerical solution and an evaluation of the parallel aspects will follow in later work.