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.