  We study $q$-variation inequality for bilinear averaging operators over
convex bodies $(G_t)_{t>0}$ defined by \begin{align*}
\mathbf{A}_t^G(f_1,f_2)(x) & =\frac{1}{|G_t|}\int_{G_t} f_1(x+y_1)f_2(x+y_2)\,
dy_1\, dy_2, \quad x\in \Bbb R^d. \end{align*} where $G_t$ are the dilates of a
convex body $G$ in $\Bbb R^{2d}$. We prove that
$$\|V_q(\mathbf{A}_t^G(f_1,f_2): t>0) \|_{L^p} \lesssim \|f_1\|_{L^{p_1}}
\|f_2\|_{L^{p_2}},$$ for $2<q<\infty$, $1<p_1,p_2\le \infty$, $1/2<p<\infty$
with $1/p=1/p_1+1/p_2$. The target space $L^p$ should be replaced by
$L^{p,\infty}$ for $p_1=1$ and/or $p_2=1$, and by dyadic BMO when
$p_1=p_2=\infty$. As applications, we obtain variational inequalities for
bilinear discrete averaging operators, bilinear averaging operators of
Demeter-Tao-Thiele, and ergodic bilinear averaging operators. As a byproduct,
we also obtain the same mapping properties for a new class of bilinear square
functions involving conditional expectation, which are of independent interest.
