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, you can pick $\omega$ to satisfy the gauge condition $d\omega=0$, and so demanding $du=F_{0}$ implies that $d\delta\omega=\Delta \omega=F_{0}$, so $\omega$ is completely determined by the constraint $du=F_{0}$. $\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$ as in my original answer.