The linearization of the operator $Pu(x)=\mathrm{div}(u)(x)-f(x,u(x))$ around $u_{0}$ is the linear operator $(D_{u_{0}}P)v=\mathrm{div}(v)(x)-df(x,u_{0}(x))v$. This linear operator has a surjective symbol, and so the operator $P$ is what is called an underdetermined elliptic operator. If $u_0$ is a function such that at $x_0$ it holds that $Pu_{0}(x_0)=0$, then $u_0$ is called an infinitesimal solution of the equation $Pu=0$ at $x_0$, and there is a local solution for $Pu=0$ around $x_0$ that asymptotically coincides with $u_0$. For reference on this, see for example Chapter 14.3 in Taylor's Partial Differential Equations, Part 3.
Since the system is underdetermined elliptic, in general the local solution will be neither unique nor sufficiently regular. A way to make $P$ into an elliptic operator, hence with certain uniqueness and regularity clauses, is to do what @Willie Wong suggested in the comments and introduce the constraint $u=\nabla w$ for a scalar function $w$. This turns the equation into a second order, nonlinear elliptic equation, and asymptotically the solution ill be unique. In your case, the linearization become $(D_{w_{0}}P)q=\Delta q(x)-df(x,\nabla w_{0}(x))\nabla q(x)$, which by restricting to an appropriate subdomain $\Omega\subset\mathbb{R}^{n}$ and imposing boundary conditions, for instance Dirichlet ones, is invertible for all $w_{0}$. If you restrict $P$ to a certain Banach space, you can then use the implicit function theorem and conclude some invertibility for $P$ when restricted thus, hence uniqueness of a solution to $Pw=0$.
Regularity in this case becomes an issue, but there are regularity theorems for nonlinear elliptic equations which may help. You can read more on this, for example, in Taylor's book.
EDIT: Demanding $u$ to be a gradient is equivalent to demanding the constraint $du=0$ (since we solve locally anyway), where $d$ is the exterior derivative, given in coordinates by $(du)_{ij}=\partial_{i}u_{j}-\partial_{j}u_{i}$. More generally, imposing a condition $du=F_{0}$ for a fixed $F_{0}$ that satisfied $dF_{0}=0$ also turns the equation into an elliptic one: by the Hodge decomposition, every vector field is of the form $u=\nabla w+\delta\omega+\kappa$ where $\kappa$ is a harmonic vector field and $\omega_{ij}$ is a $2$-form such that $(\delta\omega)_{j}=-\sum_{i}\partial_{i}\omega_{ij}$ (we work in a Euclidean setting, so I move from vector fields to differential forms indiscriminately. The arguments are the same in a Riemannian setting by using the musical isomorphism).
Now, demanding $du=F_{0}$ implies that $d\delta\omega=F_{0}$. The condition $dF_{0}=0$ implies that you can solve this equation uniquely for $\omega$, so it effectively insures that $\omega$ is completely determined by the constraint $du=F_{0}$.
As for the other terms in the Hodge decomposition of $u$, $\kappa$ is a finite dimensional obstruction that vanishes if you assume that you solve the equation $Pu=0$ locally. You are left with determining the scalar function $w$. However, since $\mathrm{div}\circ\delta=0$, you return to the situation where you effectively replace $u=\nabla w$, turning the equation $Pu=0$ into an ellitpic second order equation for $w$:
$Pw(x)=\Delta w(x)-f(x,\nabla w(x)+\delta\omega(x))=0$
Note how $\omega$ is fixed, and the only variable is $w$, and so you can deal with this equation in the same manner as my original answer.
It is worth to note that since the Hodge decomposition is valid for a differential form of any order, this disucssion applies also if you replace $u$ with a differential form of any order, with appropriatly replacing $\nabla$ with the exterior derivative $d$ and $\mathrm{div}$ with the dual of $d$, denoted by $\delta$.
