Basic FEM knowledge should tell yout that Neumann boundary conditions end up being a boundary integral on the right hand side of the Neumann datum times the test function. Homogeneous Neumann ends up being a boundary integral of zero times test, hence zero. Therefore: drop the application of the Dirichlet BC, and don’t add anything on the rhs.

I’ve run the simulation without the Dirichlet boundary condition for the top wall, and the simulation is no longer stable. I imagine that I need to specify another condition since I’ve removed a constraint.

As the code you are referring to is using a splitting scheme, you need to apply the same kind of conditions on the top surface as is done on the outlet boundary (as you are more or less applying the same constraint).