Problem on convergence rate for Reynolds no=10 or above 10, for time dependent Navier-Stokes equations

Dear all,
I am trying to compute the L2 velocity rate for time dependent Navier-Stokes equations using P2-P1 finite element space. when I take Reynolds number =1, then it works well, but when I take Reynolds number= 10 or above 10, then the L2 error convergence rate for velocity decreases and it is not coming 3, according to theory. What should I do to recover the convergence rate?

Here I have attached the code for Reynolds number =10.

/*

  • Incompressible Navier Stokes with Taylor-Hood Finite element
  • NSE use Newton method for solution 1
  • Sol x y
    */
    //load “Element_P3”

verbosity = 0;
// Parameters
real nu = 1./10.;

int maxiter = 6; // Maximum iteration
real[int] L2errorV(maxiter);
real[int] H1errorV(maxiter);
real[int] L2errorP(maxiter);
real[int] cputime(maxiter); // for store cpu time

// Start mesh loop
for (int i=1; i<= maxiter; i++ ) {
real cpu = clock();
// mesh size
int NN = 2^(i); // mesh size
real dt = 1./NN^3;
real T = 1; // final time
real TN = T/dt; // Total number of time iteration
real eps = 1e-6; // newton error
real epsilon = 1e-6; // // penalty parameter

// Mesh
mesh Th = square(NN, NN, flags=0);

// Fespace
fespace Xh(Th, P2);
Xh u1, u2, du1, du2, v1, v2;
Xh up1, up2, uu1, uu2, psi, phi;

fespace Mh(Th, P1);
Mh p, q, pp, pp1, dp;

// Macro
macro div(u1, u2) (dx(u1) + dy(u2)) //
//macro grad(u1, u2) [dx(u1), dy(u2)] //
macro ugrad(u1, u2, v) (u1*dx(v) + u2*dy(v)) //
macro Ugrad(u1, u2, v1, v2) [ugrad(u1, u2, v1), ugrad(u1, u2, v2)] //

// Define exact solution as a func of x y and t 
func real u1exact(real t)
		{
			return exp(t)*2*x^2*(x-1)^2*y*(y-1)*(2*y-1);
		}
func real u2exact(real t)
		{
			return -exp(t)*2*x*(x-1)*(2*x-1)*y^2*(y-1)^2;
		}
func real pexact(real t)
		{
			return exp(t)*2*(x-y);
		}

// Define right hand side for steady Stokes
Xh fL1 =  - nu* (4*y*(2*y - 1)*(x - 1)^2*(y - 1) + 4*x^2*y*(2*y - 1)*(y - 1)
                     + 8*x*y*(2*x - 2)*(2*y - 1)*(y - 1) + 8*x^2*(x - 1)^2*(y - 1) + 4*x^2*(2*y - 1)*(x - 1)^2 + 8*x^2*y*(x - 1)^2)
					 + 2;
Xh fL2 =  - nu* (- 8*y^2*(x - 1)*(y - 1)^2 - 4*y^2*(2*x - 1)*(y - 1)^2 - 8*x*y^2*(y - 1)^2 
              - 4*x*(2*x - 1)*(x - 1)*(y - 1)^2 - 4*x*y^2*(2*x - 1)*(x - 1) - 8*x*y*(2*x - 1)*(2*y - 2)*(x - 1))
			  -2;
			  
// Right hand side for unsteady NSE
func real f1(real t)
		{
			return exp(t)*2*x^2*(x-1)^2*y*(y-1)*(2*y-1) - nu* exp(t)*(4*y*(2*y - 1)*(x - 1)^2*(y - 1) + 4*x^2*y*(2*y - 1)*(y - 1)
                     + 8*x*y*(2*x - 2)*(2*y - 1)*(y - 1) + 8*x^2*(x - 1)^2*(y - 1) + 4*x^2*(2*y - 1)*(x - 1)^2 + 8*x^2*y*(x - 1)^2)
					 + 2*exp(t) + exp(t)*2*x^2*(x-1)^2*y*(y-1)*(2*y-1) *exp(t)*4*x*(x-1)*(2*x-1)*y*(y-1)*(2*y-1) 
					 -exp(t)*2*x*(x-1)*(2*x-1)*y^2*(y-1)^2* exp(t)* 2*x^2*(x-1)^2*((2*y-1)^2+2*y*(y-1));
		}
func real f2(real t)
		{
			return -exp(t)*2*x*(x-1)*(2*x-1)*y^2*(y-1)^2 - nu* exp(t)*(- 8*y^2*(x - 1)*(y - 1)^2 - 4*y^2*(2*x - 1)*(y - 1)^2 - 8*x*y^2*(y - 1)^2 
					- 4*x*(2*x - 1)*(x - 1)*(y - 1)^2 - 4*x*y^2*(2*x - 1)*(x - 1) - 8*x*y*(2*x - 1)*(2*y - 2)*(x - 1))
					-2*exp(t) + exp(t)*2*x^2*(x-1)^2*y*(y-1)*(2*y-1) * (-exp(t)*2*y^2*(y-1)^2*((2*x-1)^2+2*x*(x-1)))
					-exp(t)*2*x*(x-1)*(2*x-1)*y^2*(y-1)^2 * (-exp(t)*4*x*(x-1)*(2*x-1)*y*(y-1)*(2*y-1));
		}


// Problem Stokes (with solve)
solve Stokes ([u1, u2, p], [v1, v2, q], solver=sparsesolver)
	= int2d(Th)
	(
	nu*( dx(u1)*dx(v1) + dy(u1)*dy(v1)
	+ dx(u2)*dx(v2) + dy(u2)*dy(v2) )
	- epsilon* p * q 				// Penalty term
	- p*div(v1, v2) - div(u1, u2)*q  )
	- int2d(Th)([fL1,fL2]'*[v1,v2])
	+ on(1, 2, 3, 4, u1=0, u2=0);
	
// Update sol for NSE
uu1 = u1;
uu2 = u2;
pp1 = p;
//plot(coef=0.2, cmm="Stokes [u1, u2] and p" ,p, [uu1, uu2], wait=1);
//cout << "Stokes L2 Error= " << sqrt(int2d(Th)((u1-u1exact)^2 + (u2-u2exact)^2)) << endl;


// Navier Stokes Solve
// Time loop

for (int j = 1; j <= TN; j++) {
real err=0;
// Update right hand side at each time step
	real t = j*dt;
	
// Newton method 
	for (int n = 0; n <= 20; n++) {

		// Solve
		solve NSE ([du1, du2, dp], [v1, v2, q], solver=sparsesolver)  
		= int2d(Th)(
		 1./dt*(du1*v1 + du2*v2)
		 +nu * (
			  dx(du1)*dx(v1) + dy(du1)*dy(v1)
			+ dx(du2)*dx(v2) + dy(du2)*dy(v2)
			)
		+ Ugrad(du1, du2, u1, u2)'*[v1, v2]
		+ Ugrad(u1, u2, du1, du2)'*[v1, v2]
		- dp*div(v1, v2) - div(du1, du2)*q
		- epsilon*dp * q 
		)
		+ int2d(Th)(
		1./dt*(u1*v1 + u2*v2)
		+ nu * (
			  dx(u1)*dx(v1) + dy(u1)*dy(v1)
			+ dx(u2)*dx(v2) + dy(u2)*dy(v2)
			  )
		+ Ugrad(u1, u2, u1, u2)'*[v1, v2]
		- p*div(v1, v2) - div(u1, u2)*q
		)
		- int2d(Th)([f1(t),f2(t)]'*[v1,v2])
		- int2d(Th) (
			1./dt*(uu1*v1 + uu2*v2)
		)
		+ on(1, 2, 3, 4, du1=0, du2=0);
		
		// Update 
		u1[] += du1[];
		u2[] += du2[];
		p[] += dp[];

		real Lu1=u1[].linfty, Lu2=u2[].linfty, Lp=p[].linfty;
		err = du1[].linfty/Lu1 + du2[].linfty/Lu2 + dp[].linfty/Lp;

		cout << "Mesh = " << NN << ", Time = " << j*dt <<  ", Newton iteration = " << n+1 << ", err = " << err << ", eps =  " << eps << endl;
		if(err < eps) break; //converge
	}
	// update NSE solution for next step
	uu1 = u1;
	uu2 = u2;
	
	// plot
	plot(coef=0.2, cmm="Time= "+t+", [u1, u2] and p" ,p, [u1, u2], wait=0);
	//plot(coef=0.2, cmm="u1" ,u1, wait=1);
	//plot(coef=0.2, cmm="u2" ,u2, wait=1);
	//plot(coef=0.2, cmm="p" ,p, wait=1);
	//plot(coef=0.2, cmm="pexact" ,pexact, p, wait=1);

} 
// Exact sol at final time
Xh u1ex = u1exact(T);
Xh u2ex = u2exact(T);
Mh pex = pexact(T);
// error 
L2errorV[i-1] = sqrt(int2d(Th)((u1-u1ex)^2 + (u2-u2ex)^2));
H1errorV[i-1] = sqrt(int2d(Th)( (dx(u1)-dx(u1ex))^2+(dy(u1)-dy(u1ex))^2 + (dx(u2)-dx(u2ex))^2+(dy(u2)-dy(u2ex))^2));
L2errorP[i-1] = sqrt(int2d(Th)((p-pex)^2));
//cout << "L2 Error= " << L2error << endl;
// cpu time
cputime[i-1] = clock()-cpu;

}

// Print error
//
cout << "L2 Error= " << L2errorV << endl;
cout << "L2 Error= " << H1errorV << endl;
cout << "L2 ErrorP= " << L2errorP << endl;

// calculte rate
real[int] rcvl2(maxiter-1);
real[int] rcvh1(maxiter-1);
real[int] rcpl2(maxiter-1);

for(int j=1; j<maxiter; j++)
{
rcvl2[j-1] = log(L2errorV[j-1]/L2errorV[j])/log(2.);
rcvh1[j-1] = log(H1errorV[j-1]/H1errorV[j])/log(2.);
rcpl2[j-1] = log(L2errorP[j-1]/L2errorP[j])/log(2.);
}

cout << "L2 Error rate vel= " << rcvl2 << endl;
cout << "H1 Error rate vel= " << rcvh1 << endl;
cout << "L2 Error rate pre= " << rcpl2 << endl;

cout << " CPU time = " << cputime << endl;