First we look for a distributional solution. Remember that, as an distribution, $\Delta u$ is defined by $$\langle\Delta u,v \rangle=-\langle\nabla u,\nabla v\rangle,\ \forall\ v\in C_0^\infty(\Omega). \tag{1}$$
From $(1)$, we can say that a solution in the distributional sense, is a function $u\in H^1(\Omega)$ with $Tu=g$ satisfying $$\int_\Omega \nabla u\nabla v=0,\ \forall\ v\in C_0^\infty(\Omega). \tag{2}$$
By density we may conclude from $(2)$ that $$\int_\Omega \nabla u\nabla v=0,\ \forall v\in H_0^1(\Omega). \tag{3}$$
This is not the only weak formulation for this problem, however, it is the one which comes from a variational problem, to wit, let $F:\{u\in H^1(\Omega):\ Tu=g \}\to \mathbb{R}$ be defined by $$Fu=\frac{1}{2}\int_\Omega |\nabla u|^2.$$
Note that $(3)$ can be rewritten as $$\langle F'(u),v\rangle=0,\ \forall\ v\in H_0^1(\Omega). \tag{4}$$
Also note that, once $\{u\in H^1(\Omega):\ Tu=g \}$ is a closed convex set of $H^1(\Omega)$ and $F$ is a coercive, weakly lower semi continuous function, we have that $F$ has a unique global minimum which satisfies $(4)$.
By not considering any kind of derivative of $u$, you can also use another weak formulation: let $C_0^{1,\Delta}(\overline{\Omega})=\{u\in C_0^1(\overline{\Omega}):\ \Delta u \in L^\infty(\Omega)\}$. A "very" weak solution, is a function $u\in L^1(\Omega)$ satisfying $$\int_\Omega u\Delta v=-\int_{\partial\Omega}g\frac{\partial v}{\partial \nu },\ \forall\ v\in C_0^{1,\Delta}(\overline{\Omega}).$$
In your setting, I mean, when $g\in H^{1/2}(\Omega)$, it can be proved that both definitions are equivalent. For references, take a look in the paper Elliptic Equations Involving Measures from Veron. It has a PDF version here. Take a look in page 8.
To conclude, I would like to adress @JLA, which gave a comment in OP; in the end, what we really want is a $H^2$ function (or more regular), because we are working with the Laplacean and it is natural to have two derivatives.
It can be proved, by using regularity theory, that $u$ is in fact in $H^2$, however, there is a huge difference between proving that $u$ is in $H^2$ after finding it in $H^1$ by the above methods and finding directly $u\in H^2(\Omega)$ by another method. Note, for example, that none of the methods above, does apply if we change $H_0^1(\Omega)$ by $H^2(\Omega)$.
After taking $v \in H^1(\Omega)$:
$$\int_{\Omega} \nabla u \cdot \nabla v dx- \int_{\partial \Omega} v (\nabla u \cdot n)ds_x = \int_{\Omega} vdx$$
Now, using the additivity of the integral on $\partial \Omega$ and the boundary conditions on the unit square:
$$\int_{\partial \Omega} v (\nabla u \cdot n)ds_x = \int_{(0,1) \times \{0\} \cup (0,1) \times \{1\}} v \partial_x u \cdot n ds_x + \int_{\{1\} \times(0,1)}v \partial_y u \cdot n_2 ds_x + \int_{\{0\} \times (0,1)} v (c(x) u(x) + \partial_y u \cdot n_2 )ds_x$$
Best Answer
First of all, I believe you have made a small sign error - the Hemholtz equation is usually given as $$\tag{$\ast$}-\Delta u -k^2u=f \qquad \text{in } \Omega. $$ (Think of eigenvalues: $-\Delta u = k^2 u$). Of course you can consider the equation as you have written it, but the you will always find that it is coercive (independent of the value of $k$) via the Poincare inequality.
The bilinear form corresponding to $(\ast)$ is: $$B(u,v) = \int_\Omega \nabla u \cdot \nabla v \, dx - k^2\int_\Omega uv \, dx. $$ The proof that $B$ is continuous will be identical to what you already have. For coercivity, \begin{align*} B(u,u) = \int_\Omega \vert \nabla u \vert^2 \, dx - k^2 \int_\Omega u^2 \, dx. \end{align*} It is well known that the first Dirichlet eigenvalue of $-\Delta$ (which is by definition the smallest) can be written as $$\lambda_1(\Omega) = \inf_{v\in H^1_0(\Omega), \, v\neq 0} \frac{\int_\Omega \vert \nabla v \vert ^2 \, dx}{\int_\Omega v^2 \, dx}.$$ Hence, for all $u\in H^1_0(\Omega)$, $$\int_\Omega u^2 \, dx\leqslant \frac 1 {\lambda_1(\Omega)} \int_\Omega \vert \nabla u \vert ^2 \, dx$$ and so $$ B(u,u) \geqslant \bigg (1-\frac{k^2}{\lambda_1(\Omega)}\bigg )\int_\Omega \vert \nabla u\vert^2 \, dx. $$ Finally, to get the $\|\cdot\|_{H^1(\Omega)}$ norm on the right hand side, we use that $\|\nabla v\|_{L^2(\Omega)}\geqslant C \|v\|_{H^1(\Omega)}$ for all $v\in H^1_0(\Omega)$ which is true via the Poincare inequality. Thus, $$B(u,u) \geqslant C\bigg (1-\frac{k^2}{\lambda_1(\Omega)}\bigg )\|v\|_{H^1(\Omega)}, $$ so $B$ is coercive provided $\lambda_1(\Omega)>k^2$ as required.