Thanks, Shaho. I had another thought about the problem and here is what I come up with for 1D line elements in a 2D space using basic trigonometry. It is quite a bit different from your code and does not solve the issue that, possible, the calculation of the interpolation function derivatives for 1D line elements oriented in 2D space may be wrong. At least if you run my simple test above where you get different pressure gradients than the true solution.

Consider the following situation

Where the red vector is the velocity within a 1D fracture that has the same orientation. The fracture is represented as a 1D line element with starting point $P1$ (at location x1, y1 and hydrostatic pressure $\pi1$) and endpoint $P2$ (at location x2, y2 with hydrostatic pressure $\pi2$). For simplicity assume that permeability and viscosity are one and only a single phase fluid is present, which is at perfect hydrostatic equilibrium such that no fluid moves up or down and $\partial p/\partial z = \rho g$.

Now, because the fracture is oriented at some angle $\alpha$ to the vertical gravity vector $\mathbf{g}$, there will be a pressure gradient in the x and z direction, causing an "apparent velocity" $\mathbf{v}$ with an x and z component. Obviously, the z component $vz$ vanishes if we subtract $\rho g$ but the x component $vx$ remains, causing an apparent flow to the right, which clearly is non physical.

If we know $\alpha$ , we can compute vector $\mathbf{x}$ using $\tan\alpha=\mathbf{x}/\mathbf{g}$. Note that $\mathbf{x}$ is equivalent to the apparent $vx$ that we would like to subtract from $\mathbf{v}$ to correct for all hydrostatic effects, both in z and x direction. Any non-hydrostatic pressure effects will not be affected. We can obtain $\alpha$ from the scalar product of $\mathbf{g}$ and $\mathbf{v}$, where the orientation of $\mathbf{v}$ is given by $P1$ - $P1$. The scalar product allows us to compute the angle between two vectors as

(1)
\begin{align} \mathbf{v}\cdot\mathbf{g}=|v|\cdot|g|\cos\alpha \end{align}

Note that the scalar product only considers the length and orientation of the vectors intrinsically and can easily be solved for $\cos\alpha$ and we can compute $\alpha$ from the secant, i.e the inverse of $\cos\alpha$. Once we know $\alpha$ , we obtain $\mathbf{x}$ from $\tan\alpha=\mathbf{x}/\mathbf{g}$. Subtracting this value and $\mathbf{g}$ from the real velocity vector then corrects for gravity, not only in the vertical direction but also in the horizontal one. Some simple calculations on a piece and paper show that this corrects properly for gravity but leaves any other pressure gradients in order.

The same technique can be applied in 3D but we need another correction for the apparent offset in the y-direction.

Does this make sense?