Laplace Equation with Tangential Derivative Prescribed on the Boundary Consider the following Laplace boundary value problem (BVP)
$$\matrix{
   {{\nabla ^2}\Phi (x,y) = 0,} \hfill & { - a \le x \le a} \hfill & { - b \le y \le b} \hfill  \cr 
   {{{\partial \Phi } \over {\partial y}}(a,y) = f(y)} \hfill & {} \hfill & {} \hfill  \cr 
   {{{\partial \Phi } \over {\partial y}}( - a,y) = f(y)} \hfill & {} \hfill & {} \hfill  \cr 
   {{{\partial \Phi } \over {\partial x}}(x,b) = 0} \hfill & {} \hfill & {} \hfill  \cr 
   {{{\partial \Phi } \over {\partial x}}(x, - b) = 0} \hfill & {} \hfill & {} \hfill  \cr 
 } $$
where $f(-y)=-f(y)$. We have prescribed the tangential derivatives of $\Phi(x,y)$ on the boundary instead of the function value. 
The more general form of the problem in 2D or 3D can be written as
$$\matrix{
   {{\nabla ^2}\Phi ({\bf{x}}) = 0,} \hfill & {{\bf{x}} \in \Omega } \hfill  \cr 
   {{\nabla _S}\Phi ({\bf{x}}) = \left[ {\left( {{\bf{I}} - {\bf{n}} \otimes {\bf{n}}} \right) \cdot \nabla \Phi } \right]({\bf{x}}) = {\bf{f}}({\bf{x}})} \hfill & {{\bf{x}} \in \partial \Omega } \hfill  \cr 
 } $$
where $\Omega$ is the domain of interest, ${\partial \Omega }$ is its boundary and $\bf{n}$ is the unit vector normal to the boundary. Here, $\nabla_S\Phi$, called the surface gradient, denotes component of $\nabla\Phi$ tangential to $\partial\Omega$ and $\bf{f}$ is a vector field on $\partial\Omega$, tangent to the boundary at all points. $\bf{I}$ is the identity mapping, $\cdot$ is the scalar product and $\otimes$ is the tensor product.

Questions
1) Is the solution to this BVP unique? If NO, what is the degree of non-uniqueness?
2) Is there a relation between the solution to this BVP and the one with Dirichlet boundary conditions, i.e., when we determine the function value on the boundary?

My Thought
I don't think that the solution is unique so I was thinking to relate this in some manners to the solution of the Laplace BVP with Dirichlet boundary conditions where the function value is prescribed over the boundary since this BVP has a unique solution.
 A: I will deal only with the simplest version (rectangle) of the question.
I will not use the symetry property of $f$.
$\bf\text{I. Connection to the Dirichlet problem.}$
Suppose that $\Phi$ is a solution of the system, and let $\Psi=\frac{\partial^2\Phi}{\partial x^2} = -\frac{\partial^2\Phi}{\partial y^2}$. Since $\Phi$ is harmonic, its partial derivatives are harmonic, too; and in particular, $\Psi$ is harmonic.
Along the horizontal sides of the rectangle we have
$$ \Psi(x,\pm b) 
= \frac{\partial^2\Phi}{\partial x^2}(x,\pm b) 
= \frac{\partial}{\partial x} \left(
 \frac{\partial\Phi}{\partial x}(x,\pm b) \right)
= \frac{\partial}{\partial x}0
= 0.$$
Along the vertical sides we have
$$ \Psi(\pm a,y) 
= -\frac{\partial^2\Phi}{\partial y^2}(\pm a,y) 
= -\frac{\partial}{\partial y}\left(
= -\frac{\partial}{\partial y}(\pm a,y) \right)
= -f'(y). $$
Hence, $\Psi$ is known along the boundary. For computing $\Psi$ we have to solve the Dirichlet problem.
$\bf\text{II. Uniqueness.}$
Clearly the solution is not unique: if $\Phi$ is a solution then $\Phi+C$ also satisfies the conditions. But this is the only freedom.
We already know that $\Psi=\frac{\partial^2\Phi}{\partial x^2}$ is unique. The function
$-\frac{\partial^2\Phi}{\partial x\partial y}$ is the harmonic conjugate of $\Psi$, so $\frac{\partial^2\Phi}{\partial x\partial y}$ is determined up to a constant. Due to $\int_{-a}^a \frac{\partial^2\Phi}{\partial x\partial y} dx = \frac{\partial\Phi}{\partial y}(a,y)-\frac{\partial\Phi}{\partial y}(-a,y)=0$, that constant is determined; we have a unique $-\frac{\partial^2\Phi}{\partial x\partial y}$.
Then 
$$ \frac{\partial\Phi}{\partial x}(u,v)=f(v)+\int_{-a}^u \frac{\partial^2\Phi}{\partial x\partial y}(t,v)dt
\quad\text{and}\quad
\frac{\partial\Phi}{\partial y}(u,v)=\int_{-b}^v \frac{\partial^2\Phi}{\partial x\partial y}(u,t)dt $$
are uniquely determined.
Hence, the partial derivatives of $\Phi$ are determined, so $\Phi$ is uniquely determined up to a constant term.
$\bf\text{III. Existence.}$
If $f$ is nice enough (say $f$ is continously differentiable and $f'(\pm b)=0$) then the Dirichet problem has a solution for $\Psi$. From that, we want to re-construct the function $\Phi$.
The function $\Psi$ is the real part of some holomorphic function $g$, $\Psi(x,y)=\Re g(x+yi)$. 
Then $g$ has a holomorphic antiderivative $g_1$,
and $g_1$ also has a holomorphic antiderivative $g_2$ such that
$g_1(-a-bi)=g_2(-a-bi)=0$.
Now take two real constants $A,B$ and let
$$ \Phi(x,y) = \Re g_2(x+yi) +A(y+a)+B(x+a)(y+b) .$$
We show that the numbers $A,B$ can be chosen in such a way that $\Phi$ is a solution of the system.
Clearly, $\Phi$ is harmonic and $\Psi = \frac{\partial^2\Phi}{\partial x^2} = -\frac{\partial^2\Phi}{\partial y^2}$.
Along each horizontal side of the rectangle, $\frac{\partial\Phi}{\partial x}$ is constant, because
$$
\frac{\partial^2\Phi}{\partial x^2}(x,\pm b) = \Psi(x,\pm b) = 0.
$$
Along the bottom side this value is
$$
\frac{\partial\Phi}{\partial x}=
\frac{\partial\Phi}{\partial x}(-a,-b) = \Re g_1(-a-bi)=0.
$$
Along the top side we have
$$
\frac{\partial\Phi}{\partial x}=
\frac{\partial\Phi}{\partial x}(a,b) = \Re g_1(a+bi)+2bB.
$$
Now set $B=\frac{-\Re g_1(a+bi)}{2b}$; then we establish $\frac{\partial\Phi}{\partial x}=0$ along the top side as well. Notice that $\Phi$ is constant both along the top and the bottom sides of the rectangle.
Along each vertical side of the rectangle,
$\frac{\partial\Phi}{\partial y}-f(y)$ is constant, because
$$
\frac\partial{\partial y}\left(\frac{\partial\Phi}{\partial y}-f(y)\right) 
= \frac{\partial^2\Phi}{\partial y^2}-f'(y) = -\Psi-f'(y) = 0.
$$
Along the left-hand side this constant is
$$
\frac{\partial\Phi}{\partial y}-f=
\frac{\partial\Phi}{\partial y}(-a,-b)-f(-b)=A-f(-b).
$$
Chosing $A=f(-b)$ we get $\frac{\partial\Phi}{\partial y}=f$ along the left-hand side.
Along the right-hand side we have $\frac{\partial\Phi}{\partial y}-f=K$ with another constant $K$. By integrating along the vertical sides we can find that $K$ must be $0$:
$$
0 = 
\int_{-b}^b 0 dy =
\int_{-b}^b \left(\frac{\partial\Phi}{\partial y}(-a,y)-f(y)\right) dy =
\Phi(-a,b)-\Phi(-a,-b) - \int_{-b}^b f = \\ =
\Phi(a,b)-\Phi(a,-b) - \int_{-b}^b f =
\int_{-b}^b \left(\frac{\partial\Phi}{\partial y}(a,y)-f(y)\right) dy =
\int_{-b}^b K dy =
2bK
$$
$$
K=0.
$$
Therefore, $\Phi$ is harmonic, $\frac{\partial\Phi}{\partial x}(x,\pm b)=0$ and $\frac{\partial\Phi}{\partial y}(\pm a,y)=f(y)$, so $\Phi$
satisfies the conditions.
A: Let $g$ be the restriction of $\Phi$ to the boundary of the square.
Your boundary data determines $g$ up to a constant.
You can then take this $g$ and solve the corresponding Dirichlet problem; standard theory tells that solutions exist uniquely.
Therefore you can conclude that solutions are unique up to shifting by constants.
(The fact that shifting a solution by a constant gives another solution is immediate, since the BVP only sees derivatives of $\Phi$.)
I assume that you want your function $\Phi$ to be continuous up to the boundary.
You did not specify regularity assumptions, so I assume you work "classically enough" for this to make sense without any Sobolev spaces or such.
In general, knowing the gradient determines a function up to a constant, and now you know the boundary values up to a gradient.
(The boundary is a manifold, not a Euclidean domain, so this can get a bit technical, but the same idea is true.)
Let us see how to find $g$ in your specific case.
First, $\partial_x\Phi(x,b)=0$, so $\Phi(x,b)=C$ for all $x$ for some constant $C$.
Then, $\partial_y\Phi(a,y)=f(y)$, so
$$
\Phi(a,y)
=
\Phi(a,b)
+
\int_b^yf(t)dt.
$$
Since $g$ is continuous at the corner point $(a,b)$, you have $\Phi(a,b)=C$.
You assumed $f$ to be antisymmetric, so $\int_{-b}^bf=0$, and so $\Phi(a,-b)=C$.
Now $\partial_x\Phi(x,-b)=0$, so in fact $\Phi(x,-b)=C$ for all $x$.
And then you can use the same idea to find $\Phi$ on the remaining edge at $x=-a$.
You will get
$$
\Phi(x,\pm b)=C\quad\text{and}\\
\Phi(\pm a,y)=C+\int_b^yf(t)dt.
$$
These are your boundary values for a standard Dirichlet problem.
Alternatively, you could solve the ODE at each boundary part separately, obtaining four constants of integration.
You will then get four equations to tie these together by continuity at the corners.
One of these equations will be redundant — unless the boundary conditions are inconsistent and there is no solution.
Remark about your general formulation:
In 2D it is fine, but in 3D not.
There is no preferred tangent unit vector field.
You either need to specify which vector field you choose or look at the whole tangential gradient instead.
On the unit sphere $S^2\subset\mathbb R^3$ there is no globally defined unit tangent vector!

Added remark about the general case:
Consider the boundary value problem for a harmonic function $\Phi$ in a bounded smooth domain $\Omega\subset\mathbb R^n$, $n\geq2$, with prescribed tangential gradient on $\partial\Omega$.
Assume that the boundary $\partial\Omega$ is connected; this is true for balls, for example.
A solution need not exist; the boundary condition can be self-contradictory.
To see why, consider the problem of finding a function $f:\mathbb R^n\to\mathbb R$ whose gradient is a prescribed vector field.
No such function $f$ exists if the vector field does not satisfy some compatibility conditions.
(For $n=3$ a sufficient and necessary compatibility condition is that the curl of the vector field vanish. The situation is a bit more complicated on manifolds than on Euclidean spaces, but the need for compatibility conditions remains. There can be non-local conditions, and the keyword is cohomology.)
Suppose that the boundary conditions are indeed such that there is at least one solution.
How unique is this solution, then?
Suppose you had two solutions $\Phi$ and $\Psi$ to your BVP.
Then $u=\Phi-\Psi$ solves the same (linear!) equation and has zero tangential gradient at the boundary.
Since $\partial\Omega$ is connected, you can join any two points on it by a smooth curve on $\partial\Omega$.
You can express the difference between the values of $u$ at these endpoints in terms of an integral of the gradient along the curve.
This leads to the observation that $u$ is constant on the boundary.
Now we can solve $u$ from a Dirichlet problem with constant boundary values, and the unique solution is well known to be constant.
Since $u$ is constant throughout $\Omega$, the solutions $\Phi$ and $\Psi$ differ by a constant.
That is, solutions are unique up to additive constants also in the general case.
If $\partial\Omega$ is not connected, the situations is more complicated.
Then the same argument shows that the difference of any two solutions is constant on each component of $\partial\Omega$.
If there are $m$ components in $\partial\Omega$, then the (affine) space of solutions to your BVP has dimension $m$.
A: Using the fundamental theorem of calculus and the Green's function $G$ for $\nabla$ that is $G(r) = r/[4\pi|r|^3]$, we can see in 3d that the following must be true:
$$\begin{align*}\nabla \Phi(r) &= -\int_V (G \circ s)(r') \times [\nabla \times (\nabla \Phi)|_{r'} ]\, dV' + \int_V (G \circ s)(r') (\nabla \cdot \nabla \Phi)|_{r'} \, dV\\
&\quad - \int_{\partial V} (G \circ s)(r') \times [\hat n(r') \times\nabla \Phi|_{r'}] \, dS' + \int_{\partial V} (G \circ s)(r') [\hat n(r')\cdot \nabla \Phi|_{r'}] \, dS'\end{align*}$$
(n.b. may be slightly wrong on the signs of these integrals, but the basic forms are correct)
where $s(r') = r - r'$ and $\hat n$ is the unit normal to the boundary. According to the  identity
$\nabla \Phi  = \left( {\nabla \Phi .\hat n} \right)\hat n - \hat n \times \left( {\hat n \times \nabla \Phi } \right)$,
specifying the tangential derivative is to specify $\hat n \times (\hat n \times \nabla \Phi)$. We know that $\nabla \times \nabla \Phi = 0$ identically, and $\nabla \cdot \nabla \Phi = 0$ as specified by the problem.
That eliminates both the volume integrals, but the surface integral with the normal derivative--corresponding to Neumann boundary conditions--still remains.  So at this point, we can already conclude that the freedom to choose the normal derivative corresponds to freedom in $\nabla \Phi$.
So with the freedom to choose the value of $\Phi$ at one particular point (which we still have, as a constant of integration), the freedom to choose the normal derivative still allows for several distinct solutions (imagine selecting the value of $\Phi$ at a particular point on the boundary, and then integrating inward, normal to the boundary, to find $\Phi$ there; different values of the normal derivative will yield different values of $\Phi$).
In short, if you allow different values of the normal derivative, then the solution necessarily is not unique (aside from the ability to add or subtract a constant function).
