I have a suggestion how to set the right boundary condition (at least in Comsol). The point is that the advection diffusion equation can be written as $$ \nabla \cdot J = 0, $$ with $$ J=v\,c - D\nabla c $$ where $v$ is the velocity field and $D$ is the diffusion coefficient.

The bulk lagrangian density (weak form) is $$ \mathcal{L}_\mathrm{bulk}=-\nabla c^\dagger\cdot J $$ and the boundary contribution is $$ \mathcal{L}_\mathrm{boundary}=c^\dagger J\cdot\hat{n} $$ this is the weak form contribution to use in place of the outlet boundary condition. This avoid the ripples for the concentration field at the outlet. Here is the solution with a computed velocity field for $Pe=100$ (almost advection dominated problem)

I have a suggestion how to set the correct outlet boundary condition (I think this is corrupting the entire solution... at least in Comsol). The point is that the advection diffusion equation can be written as $$ \nabla \cdot J = 0, $$ with $$ J=v\,c - D\nabla c $$ where $v$ is the velocity field and $D$ is the diffusion coefficient.

The bulk lagrangian density (weak form) is $$ \mathcal{L}_\mathrm{bulk}=-\nabla c^\dagger\cdot J $$ where $c^\dagger$ is the test functtion, and the boundary contribution is $$ \mathcal{L}_\mathrm{boundary}=c^\dagger J\cdot\hat{n} $$ this is the weak form contribution to use in place of the outlet boundary condition. This avoid the ripples for the concentration field at the outlet. Here is the solution with a computed velocity field for $Pe=100$ (almost advection dominated problem)