## 1 Introduction  
Numerical simulation of incompressible flows governed by the incompressible Navier-Stokes equations has to deal with unique issues not present in compressible equations because of the lack of an independent equation for the pressure, whose gradient contributes to each of the three momentum equations and the incompressibility constraint expressed by the continuity equation stating that the velocity vector field is divergence-free. One of the common approaches to overcome the difficulties is by means of so-called projection methods originating from the ground-breaking work of [PERSON] and [PERSON] [1, 2]. The key idea of projection methods is that the momentum equations are time-advanced without satisfying the incompressibility constraint, then a correction is applied to the provisional velocity field to project it onto a divergence-free space. Practically, at each time step, one only needs to solve a sequence of decoupled elliptic equations for the velocity and the pressure. This attractive feature makes projection methods very efficient for large-scale numerical simulations. Projection methods can be classified into three classes, namely the pressure-correction schemes, the velocity-correction schemes, and the consistent splitting schemes [3]. Among them, pressure-correction schemes are the most-widely used. Pressure-correction schemes consist of two sub-steps for each time step: the first sub-step advances the momentum equations in time with the pressure treated explicitly or ignored, and the second sub-step performs the pressure correction in which the pressure is obtained by solving the derived pressure Poisson equation; the intermediate velocity field resulting from the first sub-step is then corrected by applying the pressure gradient so that the incompressibility constraint is enforced. Historically, both implicit (e.g., [4; 5; 6]) and explicit (e.g., [7; 8; 9; 10]) time-integration schemes have been used in the first sub-step. Implicit methods are desirable for steady-state or slow-transient flows where the use of large time steps is preferred, but require the solution of large non-linear algebraic equations besides the Poisson equation for the pressure. Explicit methods, on the other hand, have a more stringent stability restriction on the time step compared to implicit methods, but are much easier to implement and faster than implicit ones for fast-transient flows for which the time step needed to accurately capture physical phenomena is comparable to that the stability restriction allows. From a high-performance computing perspective, explicit methods are more attractive because they do not require global non-linear solvers, hence have increased arithmetic intensity (namely, favour computation over data communication and transfers). Due to these merits, pressure-correction explicit projection methods have been successfully used in many academic and industrial applications including direct numerical simulation (DNS) of laminar and turbulent flows (e.g., [11; 12; 13]), large-eddy simulation (LES) of atmospheric boundary layers (e.g., [7; 14; 15]), and numerical simulation of flow problems by meshless methods such as smoothed particle hydrodynamics (SPH) (e.g., [16; 17; 18]). Although projection methods and some other popular approaches have been well established for several decades, the quest for cost-effective numerical algorithms and efficient software implementations is continuing, especially in conjunction with the need for high-fidelity simulations using high-performance computing facilities. For pressure-correction explicit projection methods, the bottleneck of computational efficiency lies in the numerical solution of the pressure Poisson equation at each time level, which requires a global linear solver and is usually the most time-consuming part in the whole simulation pipeline. For multi-stage high-order Runge-Kutta (RK)-based projection methods, the pressure-related cost is further increased, because a Poisson solver is typically required at each intermediate RK stage.