Agree with @bro for the discontinuity issue. However, it might be much clearer if one writes down the governing equation for $u$ when $J[u]$ is minimized, instead of simply going for the expression of the first variation. In fact, the discontinuity of $\chi$ on $\Omega$ would lead to a jump condition for $\nabla u$.
Notations
Suppose the original functional
$$
J[u]=\int_{\Omega}\left(\left\|\nabla u\right\|^2+Q^2\cdot\mathbb{1}_{\left\{u>0\right\}}\right){\rm d}V
$$
is equipped with a boundary condition
$$
u|_{\partial\Omega}=g,
$$
where $\Omega\subseteq\mathbb{R}^n$, and ${\rm d}V={\rm d}x_1{\rm d}x_2\cdots{\rm d}x_n$. Suppose, in a more general case, that $Q$ is a function of $\mathbf{x}\in\Omega$ and $u$, i.e., $Q=Q(\mathbf{x},u)$.
Now, consider a family of feasible functions
$$
v:\Omega\times\mathbb{R}\to\mathbb{R},\quad\left(\mathbf{x},t\right)\mapsto v(\mathbf{x},t)
$$
with
$$
v|_{\partial\Omega}=g,
$$
where $u=v(\mathbf{x},0)$ gives the minima of the functional $J[v]$, i.e.,
$$
\left(\frac{\rm d}{{\rm d}t}J[v(\cdot,t)]\right)\Bigg|_{t=0}=0.
$$
Here $t$ is an auxiliary parameter, somehow playing the role of "time": as $t$ changes, the set
$$
A_t:=\Omega\cap\left\{v(\cdot,t)>0\right\}
$$
changes as well, making its boundary $\partial A_t$ deform continuously.
Finally, provided that $\nabla v(\cdot,0)$ may not be well defined on $\partial A_0$, let us define
\begin{align}
w(\cdot,t)&=v(\cdot,t)|_{A_t},\\
\omega(\cdot,t)&=v(\cdot,t)|_{B_t}
\end{align}
for further clarification, where $B_t=\Omega\setminus A_t$.
First variation: Formulism
With the notations above, our target functional reads
$$
J[v]=\int_{B_t}\left\|\nabla\omega\right\|^2{\rm d}V+\int_{A_t}\left(\left\|\nabla w\right\|^2+Q^2(\cdot,w)\right){\rm d}V.
$$
The first variation of $J$ is no more than
$$
\frac{\rm d}{{\rm d}t}J[v(\cdot,t)].
$$
Note that in $J$, not only functions $\omega$ and $w$ depend on $t$, but also integral domains $A_t$ and $B_t$. Thus to take the derivative with respect to $t$, one must use Leibniz integral rule. With this rule, we have
\begin{align}
\frac{\rm d}{{\rm d}t}\int_{B_t}\left\|\nabla\omega\right\|^2{\rm d}V&=\int_{B_t}\frac{\partial}{\partial t}\left\|\nabla\omega\right\|^2{\rm d}V+\int_{\partial B_t}\left\|\nabla\omega\right\|^2\mathbf{w}\cdot{\rm d}\mathbf{S}\\
&=2\int_{B_t}\nabla\omega\cdot\nabla\dot{\omega}{\rm d}V+\int_{\partial B_t}\left\|\nabla\omega\right\|^2\mathbf{w}\cdot{\rm d}\mathbf{S}\\
&=2\int_{B_t}\left[\nabla\cdot\left(\dot{\omega}\nabla\omega\right)-\dot{\omega}\Delta w\right]{\rm d}V+\int_{\partial B_t}\left\|\nabla\omega\right\|^2\mathbf{w}\cdot{\rm d}\mathbf{S}\\
&=-2\int_{B_t}\dot{\omega}\Delta\omega{\rm d}V+2\int_{\partial B_t}\dot{\omega}\nabla\omega\cdot{\rm d}\mathbf{S}+\int_{\partial B_t}\left\|\nabla\omega\right\|^2\mathbf{w}\cdot{\rm d}\mathbf{S},
\end{align}
where $\dot{\omega}$ stands for $\partial\omega/\partial t$, while $\mathbf{w}$ denotes the normal "velocity" of the moving part of $\partial B_t$ with respect to "time" $t$. We will come back to this $\mathbf{w}$ in the next section.
Similarly, we have
\begin{align}
&\frac{\rm d}{{\rm d}t}\int_{A_t}\left(\left\|\nabla w\right\|^2+Q^2(\cdot,w)\right){\rm d}V\\
&=\int_{A_t}\frac{\partial}{\partial t}\left(\left\|\nabla w\right\|^2+Q^2(\cdot,w)\right){\rm d}V+\int_{\partial A_t}\left(\left\|\nabla w\right\|^2+Q^2(\cdot,w)\right)\mathbf{w}\cdot{\rm d}\mathbf{S}\\
&=2\int_{A_t}\left(\nabla w\cdot\nabla\dot{w}+Q\frac{\partial Q}{\partial w}\dot{w}\right){\rm d}V+\int_{\partial A_t}\left(\left\|\nabla w\right\|^2+Q^2(\cdot,w)\right)\mathbf{w}\cdot{\rm d}\mathbf{S}\\
&=2\int_{A_t}\left(-\Delta w+Q\frac{\partial Q}{\partial w}\right)\dot{w}{\rm d}V+2\int_{\partial A_t}\dot{w}\nabla w\cdot{\rm d}\mathbf{S}+\int_{\partial A_t}\left(\left\|\nabla w\right\|^2+Q^2(\cdot,w)\right)\mathbf{w}\cdot{\rm d}\mathbf{S}.
\end{align}
Depiction for $\mathbf{w}$
To figure out a depiction for the normal velocity of $\partial A_t$ (which is also the moving part of $\partial B_t$), $\mathbf{w}$, consider a particle $\mathbf{y}=\mathbf{y}(t)$ that tracks $\left\{w=0\right\}$, i.e., $\mathbf{y}(t)$ is such that
$$
w(\mathbf{y}(t),t)=0.
$$
This equality leads to
$$
0=\frac{\rm d}{{\rm d}t}w(\mathbf{y}(t),t)=\dot{\mathbf{y}}(t)\cdot\nabla w(\mathbf{y}(t),t)+\dot{w}(\mathbf{y}(t),t).
$$
Suppose, in addition, that $\mathbf{y}(t)$ moves perpendicular to $\partial A_t$, for which $\dot{\mathbf{y}}(t)=\mathbf{w}(t)$. Using this assumption, the last equality reduces to
$$
\dot{w}(\mathbf{y}(t),t)+\mathbf{w}(t)\cdot\nabla w(\mathbf{y}(t),t)=0.
$$
Provided the arbitrariness of $\mathbf{y}(t)\in\partial A_t$, one eventually obtains
$$
\dot{w}+\mathbf{w}\cdot\nabla w=0
$$
holds on $\partial A_t$. Similarly,
$$
\dot{\omega}+\mathbf{w}\cdot\nabla\omega=0
$$
holds on the moving part of $\partial B_t$.
Surface integral of vector fields: Simplification
Let $\mathbf{n}$ be the outward unit normal of $A_t$.
Note that
\begin{align}
A_t&=\left\{w>0\right\},\\
\partial A_t&=\left\{w=0\right\}.
\end{align}
Thus
$$
\mathbf{n}=-\frac{\nabla w}{\left\|\nabla w\right\|}.
$$
Consequently,
$$
\mathbf{w}\cdot{\rm d}\mathbf{S}_{\partial A_t}=\left(\mathbf{w}\cdot\mathbf{n}\right){\rm d}S_{\partial A_t}=-\left(\mathbf{w}\cdot\frac{\nabla w}{\left\|\nabla w\right\|}\right){\rm d}S_{\partial A_t}=\frac{\dot{w}}{\left\|\nabla w\right\|}{\rm d}S_{\partial A_t}.
$$
Similarly,
$$
\mathbf{n}=-\frac{\nabla\omega}{\left\|\nabla\omega\right\|}.
$$
Yet since $\mathbf{n}$ is the outward unit normal of $A_t$, it is the inward unit normal of $B_t$. Consequently,
$$
\mathbf{w}\cdot{\rm d}\mathbf{S}_{\partial B_t}=-\left(\mathbf{w}\cdot\mathbf{n}\right){\rm d}S_{\partial B_t}=\left(\mathbf{w}\cdot\frac{\nabla\omega}{\left\|\nabla\omega\right\|}\right){\rm d}S_{\partial B_t}=-\frac{\dot{\omega}}{\left\|\nabla\omega\right\|}{\rm d}S_{\partial B_t}.
$$
Likewise,
\begin{align}
\nabla w\cdot{\rm d}\mathbf{S}_{\partial A_t}&=\left(\nabla w\cdot\mathbf{n}\right){\rm d}S_{\partial A_t},\\
\nabla\omega\cdot{\rm d}\mathbf{S}_{\partial B_t}&=-\left(\nabla\omega\cdot\mathbf{n}\right){\rm d}S_{\partial B_t}.
\end{align}
First variation: Results
With all the results in the surface integral of vector fields, the first variation reads
\begin{align}
\frac{\rm d}{{\rm d}t}J[v]&=2\int_{A_t}\left(-\Delta w+Q\frac{\partial Q}{\partial w}\right)\dot{w}{\rm d}V+2\int_{B_t}\left(-\Delta\omega\right)\dot{\omega}{\rm d}V+\\
&\quad\quad\int_{\partial A_t}\left(\left\|\nabla w\right\|+\frac{Q^2}{\left\|\nabla w\right\|}+2\mathbf{n}\cdot\nabla w-\left\|\nabla\omega\right\|-2\mathbf{n}\cdot\nabla\omega\right)\dot{v}{\rm d}S.
\end{align}
Provided that
\begin{align}
\mathbf{n}\cdot\nabla w&=-\frac{\nabla w}{\left\|\nabla w\right\|}\cdot\nabla w
=-\left\|\nabla w\right\|,\\
\mathbf{n}\cdot\nabla\omega&=-\frac{\nabla\omega}{\left\|\nabla\omega\right\|}\cdot\nabla\omega=-\left\|\nabla\omega\right\|,
\end{align}
the result is equivalent to
\begin{align}
\frac{\rm d}{{\rm d}t}J[v]&=2\int_{A_t}\left(-\Delta w+Q\frac{\partial Q}{\partial w}\right)\dot{w}{\rm d}V+2\int_{B_t}\left(-\Delta\omega\right)\dot{\omega}{\rm d}V+\\
&\quad\quad\int_{\partial A_t}\mathbf{n}\cdot\left[\left(1-\frac{Q^2}{\left\|\nabla w\right\|^2}\right)\nabla w-\nabla\omega\right]\dot{v}{\rm d}S.
\end{align}
Finally, due to the arbitrariness of $\dot{v}(\cdot,0)$, the governing equations are
\begin{align}
-\Delta w+Q\frac{\partial Q}{\partial w}&=0,&\text{in }A_0,\\
-\Delta\omega&=0,&\text{in }B_0,\\
\mathbf{n}\cdot\left[\left(1-\frac{Q^2}{\left\|\nabla w\right\|^2}\right)\nabla w-\nabla\omega\right]&=0,&\text{on }\partial A_0,
\end{align}
together with the continuity of $v(\cdot,0)$, i.e.,
\begin{align}
w-\omega&=0,&\text{on }\partial A_0,
\end{align}
as well as the boundary condition for $v(\cdot,0)$, i.e.,
\begin{align}
v(\cdot,0)&=g,&\text{on }\partial\Omega.
\end{align}
As per the above governing equations, one can see that the characteristic-function $\chi$ in the functional $J$ eventually contributes to the jump condition for $\nabla v(\cdot,0)$ on the level set $\partial A_0$.
