I have been studying the 2D time-independent advection equation on the unit square $[0,1] \times [0,1]$. One such example is:
$$
\frac{\partial}{\partial x} u(x,y) + \frac{\partial}{\partial y} u(x,y) = 1 \, ,
$$ 
with the (no inflow) boundary condition $u(\cdot,0) = u(0,\cdot) = 0$.


A simple analysis (method of characteristics etc.) shows that the derivatives of the solution are discontinuous. 


My intuition tells me that this is because of the non-smooth nature of the square boundary. I'm almost certain there will be at least some papers on related results, but I cannot seem to find them. Can anyone point me in the right direction?



Advection equation regularity (2D and time independent)
The problem is that **your boundary conditions are incompatible with your equation**. 


Your choice of boundary data implies that, were $u$ to be in $C^1([0,1]\times[0,1])$ (meaning that the derivative extends continuously to the boundary), you must have $\nabla u(0,0) = 0$, since the two partials both vanish. This is clearly incompatible with the equation which would require $\partial\_x u + \partial\_y u = 1$ there. 


Philosophically then the singularity is expected to persist in the interior (unlike in the case of elliptic equations) due to the lack of smoothing. 


Another way to think about this is that for the advection equation you wrote down, the Cauchy problem w.r.t. initial data prescribed on the $x$-axis is well-posed, as is the one with respect to the $y$-axis. So when you propose the Cauchy problem w.r.t. the boundary you give, at the origin your system is **over-determined** due to it solving **two** separately well-posed IVPs at the same time. So to get solvability (in $C^1$ class) of your equation you need extra compatibility conditions at the origin. 


Note that non-smoothness of your boundary is not going to *always* cause problems. For example, had you decided to prescribe the boundary data that vanishes on $\{0\}\times[0,1]$ and equals $x$ on $[0,1]\times\{0\}$, then you will have a $C^1$ (in fact $C^\infty$) solution) to your problem. More generally, you have, for your particular system:



> 
> Consider the equation $\partial\_x u + \partial\_y u = 1$ with initial data prescribed on $\Sigma = \{0\}\times [0,\infty) \cup [0,\infty) \times\{0\}$. Suppose $f:\mathbb{R}^2\to\mathbb{R}$ is a $C^1$ function such that 
>  $$ \partial\_x f + \partial\_y f |\_{\Sigma} = 1 $$
>  Then there exists a unique $C^1$ solution $u$ such that $u = f|\_{\Sigma}$. 
> 
> 
> 


This is part of a general principle about well-posed hyperbolic PDEs. Hyperbolicity really has two components: the first is a Cauchy-Kovalevskaya-esque step that states that the "initial data" is compatible with the existence of a formal solution (think of this as a Taylor polynomial), the second step is "energy estimates" that show that this formal solution can be extended to an actual solution. [The first step is also okay usually for elliptic PDEs; for them the second step usually fails for the lack of energy estimates.] In your set-up the problem is that the first step already fails. 


