Thanks Andrew, I got an answer from my friend Enrique Nadal Soriano, and I am just copying it here:

"In FEM in general the gradient of the solution is evaluated with the

gradient of the shape functions. Usually the shape functions are

evaluated in a local space and thus its gradients. When evaluate these

gradients in the physical space (X,Y,Z) coordinates we need the Jacobian

to transformate the derivatives in the local space to the global space,

that at the end is relaed with the element size. I hope to answer your

question!"

In my oppinion the use of the Jacobian is not sufficient to give a satisfing explanation to your observation. The first displacement gradient is related to the strain, the second displacement gradient, is thus the spatial gradient of the strain. Using finite elements, also higher order finite elements one still usually only requires C0 continuity, between the elements. That means that the displacement field in the model is continuous, but the strain tensor field is not, and the strain tensor gradient field even less (in particular at the element interfaces). So although higher order elements might catch the local second order displacement gradient field within the element, as long as the second order displacment gradient field is not converged such that it is approximatly continuous at the element interfaces, this is a misleading accuracy for second derivative information. To asses second order gradient accuracy the mesh should also be converged with respect to second order gradients, and not only to the usual required (first order derivative related) stress, or strain convergence.

If the problem does not contain a source for (0th, 1st, or 2nd, order) singularity in the area of interest, choosing a fine enough mesh should however decrease the influence of the element size on a local second order gradient. If you are already on the limit of available refinement, maybe you can use some c1 continuous element formulation, but they are quite rare i think.