Chemotaxis is the ability of some bacteria to direct their movement according to the gradient of chemicals contained in their environment. In soil, some bacteria microorganisms that degrade organic carbon (SOC for short) are motile and chemotactic. This phenomenon is observed in experiments (see [1]) and on field. Nevertheless to our best knowledge no model of terrestrial carbon cycle adresses this issue. Indeed, these models are essentially compartimental corresponding naturally to systems of ordinary differential equations (e.g. Century, RothC, MOMOS) (see [2]). They are used globally to estimate soil CO2 emissions in local land management and crop optimization, among other things. Very few prototypes of spatial soil organic model have been proposed. Some of them use systems of partial differential equations: Balesdent et al. [3] combined vertical directed transport of organic carbon with a degradation phenomenon and diffusion. More recently, Goudjo et al. [4] proposed a three dimensional model for dissolved organic matter using also a system of PDEs. Other authors opted for a finite sequence of interconnected systems of ordinary differential equations each localized in a soil layer (see [5]). We previously studied the model MOMOS proposed by Pansu [6–۷], which is a nonlinear system of ordinary differential equations (see [8]) written as In these equations the unknown u models the alive microbial biomass, whereas the unknowns v and w are soil organic matters with distinct decomposition rates. In reality, the nonnegative functions ki , i ∈ {۱, ۲, ۳, ۴}, q and f depend not only on time but also on space because of the variability in soil clay content. The phenomena described by MOMOS can also be subjected to the influence of transport and sedimentation through transport and diffusion. In order to test the effect of soil heterogeneity we studied in [9] the following reactiondiffusion-advection initial problem: here Ω is a domain in R n representing the soil, Ai is a diffusion matrix and Bi a transport vector, for each i = 1, 2, 3. In [9] the boundary conditions were either of Dirichlet type (γ = ۰, βi ≡ ۱) or of NeumannRobin type (γ = ۱, βi(t, x) ≥ ۰). The right hand side term of (1.1) was where we replaced the term q(t, x)u 2 1 with q(t, x)|u1|u1 for more accuracy, since q(t, x)u1 corresponds to a kinetic coefficient that cannot be negative.