I am entertaining the following scenario:
It is well known that given any surjective submersion $p\colon M\rightarrow N$ between smooth manifolds (without boundaries), one can always find a local section $\sigma$ of $p$ passing through an arbitrary point $x\in M$, i.e. $\sigma(p(x))=x$.
Suppose now that $D\subset TM$ is a (regular, possibly non-integrable) distribution on $M$, such that the restriction of the differential $$\mathrm dp\colon D\rightarrow TN$$ is fibrewise surjective. My question is: given any $x\in M$, can we find a local section $\sigma$ of $p$ through $x$, such that the image of $\mathrm d\sigma$ lies in $D$?
Clearly, this holds whenever $D$ is integrable, since then the restriction of $p$ to any leaf is a submersion, so we can use the statement from the first paragraph to get our wanted conclusion. Moreover, note that when $\mathrm{rank}(D)=\dim (N)$, the existence of such sections implies integrability of $D$ since their images are just the integral manifolds of $D$.
Any thoughts on the case when $D$ is non-integrable?
As a simple counterexample (to the conjecture that such sections should exist), consider the distribution $D$ on $\mathbb{R}^3$ spanned by $X_1 = \partial_x - y \partial_z$ and $X_2 = \partial_y$, and the projection $\mathbb{R}^3 \to \mathbb{R}^2$ to the $xy$-plane. As you say, if $D$ had local sections it would be integrable, but it clearly is not, as $[X_1, X_2] = \partial_z$.
I think (but haven't checked carefully) that one should be able to use this construction to get counterexamples for arbitrary $\dim(D) \geq \dim(N)$. As an outline: Consider the projection $p: \mathbb{R}^n \to \mathbb{R}^2$ to the $(x_1 x_2)$-plane again, and let $D$ be the distribution spanned by $$ \begin{align*} X_1 &= \partial_1 - x_2 \partial_n \\ X_i &= \partial_i \text{ for $2 \leq i \leq n-1$}. \end{align*} $$ If we can find a local section $\sigma$ of $p$ through some point $x \in \mathbb{R}^n$ with $d\sigma$ having image in $D$, then the tangent space to $\sigma(U)$ is a 2-dimensional integrable sub-distribution $D'$ of $D$ on the restriction of $T\mathbb{R}^n$ to $\sigma(U)$. (You can think of this instead as living in the pullback of $T\mathbb{R}^n$ by $\sigma$ if you like, so that it is a distribution on an $n$-dimensional vector bundle over $U$.)
$D'$ is spanned by some vector fields $Y_1 = \sum c_{1i} X_i$ and $Y_2 = \sum c_{2i} X_i$, where $c_{ij}$ are smooth functions on $\sigma(U)$; the surjectivity of $dp$ should be equivalent to the $2 \times 2$ matrix $(c_{ij})_{1 \leq i, j \leq 2}$ being nonsingular on $U$; performing a linear transformation (on the vector bundle fibers) should bring it into upper-triangular form, so we can assume $c_{11}$ and $c_{22}$ are nonvanishing on $U$ and $c_{12} = 0$; then $[Y_1, Y_2]$ will contain a $\partial_n$ term and hence can't live in $D' \subseteq D$.