Galerkin Least-Squares stabilization for FEM solution advection (hyperbolic) equations

I am playing with Galerkin Least-Squares stabilization to solve advection diffusion problem in the context of the finite element method. This works very well for steady-state advection-diffusion equation, even in the High Péclet limit, when the problems are solved in the steady state.

The equation we solve is :
$$frac{partial alpha}{partial t} + mathbf{u} cdot nabla alpha = D nabla^2 alpha $$

with $ mathbf{u}$ a velocity vector field which is known and $alpha$ a scalar. $D$ is a diffusion coefficient.
The weak form (without boundary terms) is :
$$ int phi_i frac{partial alpha}{partial t} dOmega + int phi_i mathbf{u} nabla alpha dOmega + int D nablaphi_i cdot nablaalpha dOmega = 0$$

Adding GLS stabilization leads to the following:
$$ int phi_i frac{partial alpha}{partial t} dOmega + int phi_i mathbf{u} nabla alpha dOmega + int D nablaphi_i cdot nablaalpha dOmega + int (frac{partial alpha}{partial t} + mathbf{u} cdot nabla alpha – D nabla^2 alpha ) (tau mathbf{u} cdot nabla phi_i) dOmega = 0 $$

Where the additional term that is added effectively consists in the strong form of the equation. We define the stabilization coefficient $tau$ as:

$$tau = left( left(2 frac{|mathbf{u}|}{h} right)^2+ left(4 frac{mathbf{D}}{h^2} right)^2 right)^{-frac{1}{2}}$$

where $h$ is the size of the cell which is calculated using the diameter of a sphere of equal volume to that of the cell. Consequently $tau$ has units of time.

Solving the steady-state version of the equation :
$$mathbf{u} cdot nabla alpha = D nabla^2 alpha $$

works perfectly fine for any Péclet number. Obviously the GLS term introduces artificial diffusion, but it ensures that the value remain strictly positive and bounded. As the mesh is refined, the solution gets better and better. My issue appears in transient simulations. There, I get loss of positivity and oscillations (overshoots) in my results. It is nothing too bad (nothing like regular Galerkin would give me), but for example, the advection of a step-function without diffusion leads to overshoot of alpha and the bounded character is not preserved.

I feel like something is missing for my GLS stabilization when my equations are transient, but I have not been able to find literature that discuss this problem exactly. What am I missing? Is it because of how I define my $tau$? For information, I use BDF-type schemes (BDF1 or BDF2) for time integration.

Edit 1:
Following the (Finite Element Method : Theory, Implementation and Application by Larson) I understand that a full GLS formulation would require me to add another term that would be :

$$ int (frac{partial alpha}{partial t} + mathbf{u} cdot nabla alpha – D nabla^2 alpha ) (- D nabla^2 phi_i)$$

but this term is trivially zero for P1 elements, so I do not think it is a necessary component here. At least this is what I saw in the source below:

Edit 2

• Results obtained for steady-state advection-diffusion at Pe=1000 with a coarse mesh

• Advection of a sharp tanh function. Initial condition and after first step (note the overshoot)

• Advection of a step function. Initial condition and after first step (note the overshoot)