Higher order of the convection term approximation

when I refer to the convection term, I mean this:

$$\frac{\text{D}}{\text{D} t} \vec{v} = \frac{\partial \vec{v}}{\partial t} + \vec{v} \cdot \vec{\nabla} \, \vec{v}$$

The logic behind this is as follows (at least to my understanding): let us have a velocity field at the time $t$. What happens when we move to the time $t + \delta t$? Well, that's easy, at least up to the first order in $t$. Time itself will just move to $t + \delta t$ (duh), position of a fictious particle of the flow $\vec{v}$ will move from a general position $\vec{r}$ to: $$\vec{r} \to \vec{r} + \vec{v} (\vec{r}, t) \delta t$$

Now the velocity field changes from $\vec{v} (\vec{r}, t)$ at the time $t$ to: $$\vec{v} (\vec{r}, t) \to \vec{v} (\vec{r} + \vec{v} (\vec{r}, t) \delta t, t + \delta t)$$

So the overall change is (up to the first order in $t$): $$\vec{v} (\vec{r} + \vec{v} (\vec{r}, t) \delta t, t + \delta t) - \vec{v} (\vec{r}, t) = \frac{\partial \vec{v}}{\partial \vec{r}} \cdot \vec{v} \, \delta t + \frac{\partial \vec{v}}{\partial t} \delta t = \left( \vec{v} \cdot \vec{\nabla} \, \vec{v} + \frac{\partial \vec{v}}{\partial t} \right) \delta t$$ which is an object that is being referred to as "convective term" or "convective derivative" (if you divide by $\delta t$) and so on. My goal is to obtain higher order approximations to this - those proportional to $\delta t^2$, $\delta t^3$ etc. What is a proper way to do it? I had several ideas, unfortunately, almost every one of them yields different result.

First of all, we can expand using Taylor theorem due to arguments themselves, obtaining: $$\left( \vec{v} \cdot \vec{\nabla} \, \vec{v} + \frac{\partial \vec{v}}{\partial t} \right) \delta t + \frac{1}{2} \left( v_i v_j \partial_i \partial_j \vec{v} + \frac{\partial^2 \vec{v}}{\partial t^2} \right) \delta t^2 + \frac{1}{6} \left( v_i v_j v_k \partial_i \partial_j \partial_k \vec{v} + \frac{\partial^3 \vec{v}}{\partial t^3} \right) \delta t^3 + \cdots$$

However, what about the first (position) argument of the velocity field? There is the term $\vec{r} + \vec{v} (\vec{r}, t) \delta t$ which implies on the order itself. If expanded in this direction first, we get: $$\vec{r} \to \vec{r} + \vec{v} (\vec{r}, t) \delta t + \frac{1}{2} \underbrace{\left( \frac{\partial \vec{v}}{\partial t} + \vec{v} \cdot \vec{\nabla} \, \vec{v} \right) \delta t^2}_{\text{based on the previous}}$$

This expansion gives a different result: $$\left( \vec{v} \cdot \vec{\nabla} + \frac{\partial}{\partial t} \right) \vec{v} \, \delta t + \frac{1}{2} \left[ \left( \frac{\partial}{\partial t} + \vec{v} \cdot \vec{\nabla} \right) \vec{v} \cdot \vec{\nabla} \, \vec{v} + v_i v_j \partial_i \partial_j \vec{v} + \frac{\partial^2 \vec{v}}{\partial t^2} \right] \delta t^2 + \cdots$$ (square brackets - $\nabla$ and time derivative act only on the first velocity, but I wanted to pinpoint the pattern)

I think that the second variant is somewhat "more correct" (which is a very unfortunate statement, as there is only one correct way to do stuff in this kind of math). My argument is that if we want to go beyond the linear approximation, we have to expand really everything - not only velocity due to Taylor theorem with first order (linear) argument, but also the argument itself due to previous approximation, as I did in the second variant. This is when I'm getting a little lost - should I have included the whole convective term in the argument expansion (both time derivative and $\vec{v} \cdot \vec{\nabla} \, \vec{v}$), or only time derivative, so it would've been: $$\vec{r} \to \vec{r} + \vec{v} \, \delta t + \frac{\partial \vec{v}}{\partial t} \delta t^2$$

So my question is: what is the correct second-order term and how can we easily derive it without going mad?

P.S.: my motivation is to try and use the higher-order variant in FEM simulation of Navier-Stokes equations to see if some spurious convective-induced artefacts vanish. However, it is clear that the higher the approximation is, the more nonlinear terms appear in the equation (higher overall power of velocity). Moreover, with the $\delta t^2$ term, the resulting equations are no longer first-order linear in time.

I don't think these manipulations are correct or useful.

The Navier-Stokes equations correctly model flow behavior under conditions where the fluid can be treated as a continuum. Here the fluid velocity represents a local volume average of the velocities of the constituent molecules when the mean-free path is many orders of magnitude smaller than the macroscopic length scales of interest.

The Navier-Stokes equations are no longer applicable for rarefied gas flow, but then, the representation of the velocity as a smooth function is even less appropriate -- and it is this smoothness you are relying on to introduce these higher-order derivatives.

The reason the convection term or material derivative $\displaystyle \frac{D \mathbf{v}}{Dt}$ is appropriate as it stands, is that it is exactly what is needed to enforce conservation of momentum. The convection term arises when we equate the rate-of-change of momentum of the fluid in a material control volume $V(t)$, which is moving with the flow, to the pressure, viscous and body forces acting on that volume of fluid.

That momentum rate-of-change for an incompressible fluid is is

$$\frac{d}{dt} \int_{V(t)} \rho \mathbf{v} \, dV = \int_{V(t)} \rho \left(\frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v}\right) \, dV = \int_{V(t)} \rho \frac{D \mathbf{v}}{D t} \, dV$$.

In taking the derivative under the integral over a region with a time-dependent boundary , the appearance of the material derivative on the RHS is a consequence of the Reynolds transport theorem which is a generalization Leibniz' rule to multiple dimensions.

Since momentum must be conserved for any possible choice for the control volume, we require locally that

$$\rho \frac{D \mathbf{v}}{D t} = -\nabla p + \mu \nabla^2 \mathbf{v}.$$

Unless you disagree with Newton's laws of motion, tinkering with the term on the LHS is unnecessary. Plus the nonlinearity without this is already more than enough to handle.

Perhaps what you have in mind is how to approximate or discretize the derivatives in the Navier-Stokes equations. In that case higher-order approximations in finite differences or finite elements can lead to better accuracy or improved convergence rate in numerical solutions.

Here I am referring to approximations like

$$\left.\frac{\partial^2 v}{\partial x^2}\right|_{x_i,y_j} \approx \frac{v_{i+1,j}- 2v_{ij} + v_{i-1,j}}{\delta x^2}$$