Half a year had past and it's ripe to answer the question :-)
It is fairly easy. First, since Archimedes’ Principle is essentially about hydrostatics, then $\boldsymbol u = 0$, i.e. nothing moves. Thus your equation turns into
$$\nabla p = - \rho \boldsymbol g$$
Then you have to integrate pressure over the whole body's surface to find the force exerted on the body by the fluid:
$$\boldsymbol F = \int_S p \cdot d \boldsymbol s = \int_V \nabla p \, dV = - \rho \boldsymbol g \int_V dV = - \rho V \boldsymbol g$$
$V$ is the volume of the body, $S$ is its surface. That's all. I've used Gauss theorem here:
$$\int_V \left(\nabla\cdot\mathbf{F}\right)dV=\int_S(\mathbf{F}\cdot\mathbf{n}) \, dS$$
However as you see I've integrated $\nabla p$ inside the body, just like there was fluid and no body. If you are a mathematician, I think you may say that given $p(\boldsymbol x)$ satisfying the hydrostatic equation outside the body, there is a unique smooth continuation of $p$ inside the body, satisfying the same equation.
But I'm rather a physician and I would say, that we can substitute the body by a bubble full of fluid and the forces acting on that bubble from the outside fluid won't change. That's just because these forces depend only on the shape of the body.
UPD Oops, I've missed the minus sign, I've corrected it --- the force exerted on the body is opposite to $\boldsymbol g$.
Let us consider the initial-value problem $\rho(x,0) = \rho_0(x)$. The method of characteristics consists in introducing the parameterization $\rho(x(t),t)$ of the density. Differentiation w.r.t. time gives
$$
\frac{\text d }{\text d t}\rho = \rho_x\, \frac{\text d }{\text d t}x + \rho_t \, .
$$
Using the PDE $\rho_t + u \rho_x = 0$, one writes the system of characteristic equations
$$
\frac{\text d }{\text d t}\rho = 0\, , \qquad
\frac{\text d }{\text d t}x = u\, .
$$
Letting $\rho(0) = \rho_0(x_0)$ gives $\rho(x(t),t) = \rho_0(x_0)$, and letting $x(0) = x_0$ gives $x(t) = x_0 + \int_{0}^{t}u(s)\,\text d s$.
Combining both equations and going back to the original variables, we have
$$
\rho(x,t) = \rho_0\left(x - \int_{0}^{t}u(s)\,\text d s\right) .
$$
The constant velocity case $u(t) = \bar u$ is included in this formula, which gives $\rho(x,t) = \rho_0(x - \bar u t)$.
Best Answer
Let's linearize the problem to get to the heart of it. (This is the same as analyzing the motion of the fluid along a very small timescale and very small length scale.)
For convenience center the timescale about $t=0$ and the length scale about $\mathbf{x}=0$.
Then: $$ \rho(\mathbf{x}, t) = \rho(0, 0) + \nabla\rho\cdot (\mathbf{x} - t\mathbf{u}) $$
Why the minus sign? Because the density at position $\mathbf{x}$, time $t$ is carried forward from the density at position $\mathbf{x} - t\mathbf{u}$, time $0$.
In some sense, it may be messing with your intuition because we're used to derivatives as "final minus initial" but in this case the derivative is actually "in minus out", ie "upstream minus downstream".