  Let $L$ be a Neumann operator of the form $L=-\Delta_{N}$ acting on
$L^2(\mathbb R^n)$. Let ${BMO}_{\Delta_{N}}(\mathbb R^n)$ denote the BMO space
on $\mathbb R^n$ associated to the Neumann operator $\L$. In this article we
will show that a function $f\in { BMO}_{\Delta_{N}}(\mathbb R^n)$ is the trace
of the solution of $${\mathbb L}u=u_{t}+L u=0, u(x,0)= f(x),$$
  where $u$ satisfies a Carleson-type condition \begin{eqnarray*}
  \sup_{x_B, r_B} r_B^{-n}\int_0^{r_B^2}\int_{B(x_B, r_B)} |\nabla u(x,t)|^2
{dx dt } \leq C <\infty, \end{eqnarray*} for some constant $C>0$. Conversely,
this Carleson condition characterizes all the ${\mathbb L}$-carolic functions
whose traces belong to the space ${BMO}_{\Delta_{N}}(\mathbb R^n)$. This result
extends the analogous characterization founded by E. Fabes and U. Neri in
({Duke Math. J.} {42} (1975), 725-734) for the classical BMO space of John and
Nirenberg. Furthermore, based on the characterization of
${BMO}_{\Delta_{N}}(\mathbb R^n)$ space mentioned above, we prove global
well-posedness for parabolic equations of Navier-Stokes type with the Neumann
boundary condition under smallness condition on intial data $u_{0}\in {{
BMO}_{\Delta_{N}}^{-1}(\mathbb R^n)}$, which is motivated by the work of P.
Auscher and D. Frey ({J. Inst. Math. Jussieu} {16(5)} (2017), 947-985).
