If a separately continuous function $f : [0,1]^2 \to \mathbb{R}$ vanishes on a dense set, must it vanish on the whole set?

The answer is yes. It appears that this was first proved by Sierpiński in 1932.

W. Sierpiński, Sur une propriété de fonctions de deux variables réelles, continues par rapport à chacune de variables, Publ. Math. Univ. Belgrade 1 (1932), 125–128. http://elib.mi.sanu.ac.rs/files/journals/publ/1/8.pdf

Here's Sierpiński's proof, translated freely into English, and using your notation.

Let $E$ be dense in $[0,1]^2$. Suppose $f$ is separately continuous and vanishes on $E$. Let $(x_0, y_0) \in \mathbb{R}^2$ be arbitrary and fix $\epsilon > 0$. By separate continuity along $y=y_0$, there exists a number $\delta_0$ such that if $|x-x_0|\le \delta_0$ then $|f(x,y_0) - f(x_0, y_0)| < \epsilon$. Now recursively construct sequences $x_n, y_n, \delta_n$ as follows. Given $(x_{n-1}, y_{n-1}, \delta_{n-1})$, choose $(x_n, y_n) \in E$ such that $|x_n - x_{n-1}| < \frac{1}{2} \delta_{n-1}$ and $|y_n - y_0| < \frac{1}{n}$; this is possible by the density of $E$. By separate continuity along $y=y_n$, choose $\delta_n$ such that for all $x$ with $|x-x_n| \le \delta_n$ we have $|f(x,y_n) - f(x_n, y_n)| < \epsilon$. Without loss of generality, we can also assume $\delta_n < \frac{1}{2} \delta_{n-1}$. In particular, for any $i \ge k$, $\delta_i \le 2^{-(i-k)} \delta_k$.

Now note that for any $n > k$, we have $$|x_n-x_k| \le \sum_{i=k}^{n-1} |x_{i+1} - x_{i}| \le \sum_{i=k}^\infty \frac{1}{2}\delta_{i} \le \frac{1}{2} \delta_k \sum_{i=k}^\infty 2^{-(i-k)} = \delta_k.$$

In particular, $x_n$ is Cauchy, so it converges to some $\xi \in [0,1]$. For any $k$ we have $|x_n - x_k| \le \delta_k$ for all $n > k$, so letting $n \to \infty$ we also have $|\xi - x_k| \le \delta_k$. Therefore $|f(\xi, y_k) - f(x_k, y_k)| < \epsilon$. But $(x_k, y_k) \in E$ so $f(x_k, y_k) = 0$ and we have $|f(\xi, y_k)| < \epsilon$. As $k \to \infty$ we have $y_k \to y_0$ so by separate continuity along $x=\xi$ we have $f(\xi, y_k) \to f(\xi, y_0)$. Thus $|f(\xi, y_0)| \le \epsilon$. Since $|\xi - x_0| \le \delta_0$, we have $|f(\xi, y_0) - f(x_0, y_0)| < \epsilon$. So we have $|f(x_0, y_0)| < 2\epsilon$. Since $\epsilon$ was arbitrary we conclude $f(x_0, y_0) = 0$.

It looks to me like the proof would go through without change if we replace $[0,1]^2$ with $X \times Y$ where $X,Y$ are metric spaces and at least one of them is complete. We can also replace the codomain $\mathbb{R}$ by any metric space $Z$.

I don't know what happens if neither $X$ nor $Y$ is complete. It would also be interesting to ask what happens for more general topological spaces.