This post presents a powerful method of reasoning that avoids a great deal of algebra and calculation. To those familiar with this method, the work is so automatic and natural that one's initial response to a question like this is "it's obvious!" But maybe it's not so obvious until you have seen the method. Therefore, all the details are explained, step by step.
Background
There are several formulas for the variance of data $\mathbf{x}=x_1, x_2, \ldots, x_n$ (with mean $\bar x = (x_1+\cdots + x_n)/n$), including
$$\operatorname{Var}(\mathbf{x}) = \frac{1}{n}\sum_{i=1}^n (x_i - \bar x)^2 = \frac{1}{n}\left(\sum_{i=1}^n x_i^2\right) - \bar x^2.\tag{1}$$
This determines the covariance of paired data $(x_1,y_1), \ldots, (x_n, y_n)$ via
$$\operatorname{Cov}(\mathbf{x}, \mathbf{y}) = \frac{1}{4}\left(\operatorname{Var}(\mathbf{x}+\mathbf{y}) - \operatorname{Var}(\mathbf{x}-\mathbf{y})\right).$$
The formula implied in the referenced covariance-with-crayons post is
$$C(\mathbf{x}, \mathbf{y}) = \sum_{i=1}^{n-1} \sum_{j=i+1}^n (x_j - x_i)(y_j - y_i) = \frac{1}{2}\sum_{i,j=1}^n (x_j - x_i)(y_j - y_i).\tag{2}$$
That post asserts $C$ is proportional to the covariance. The constant of proportionality $c(n)$ could (and does) vary with $n$. Thus, when $\mathbf{x}=\mathbf{y}$ one implication of this assertion is that
$$C(\mathbf{x}, \mathbf{x}) = c(n) \operatorname{Var}(\mathbf{x}).$$
Analysis
Although this could be demonstrated with brute-force algebra, there's a better way: let's exploit the fundamental properties of covariance. Which properties would those be? I would like to suggest the following are basic:
Location invariance. That is, $$\operatorname{Cov}(\mathbf{x}, \mathbf{y}) = \operatorname{Cov}(\mathbf{x}-\mathbf{a}, \mathbf{y})$$ for any number $a$. (The expression $\mathbf{x}-\mathbf{a}$ refers to the dataset $x_1-a, x_2-a, \ldots, x_n-a$.)
Multilinearity. This implies $$\operatorname{Cov}(\lambda\,\mathbf{x}, \mathbf{y}) = \lambda\,\operatorname{Cov}(\mathbf{x}, \mathbf{y})$$ for any number $\lambda$. (The expression $\lambda\mathbf{x}$ refers to the dataset $\lambda x_1, \lambda x_2, \ldots, \lambda x_n$.)
Symmetry. The covariance of $\mathbf{x}$ and $\mathbf{y}$ is the covariance of $\mathbf{y}$ and $\mathbf{x}$: $$\operatorname{Cov}(\mathbf{x}, \mathbf{y}) =\operatorname{Cov}(\mathbf{y}, \mathbf{x}).$$
Invariance under permutations. The covariance does not change when we re-index the $(x_i, y_i)$. Formally, $$\operatorname{Cov}(\mathbf{x}, \mathbf{y}) = \operatorname{Cov}(\mathbf{x}^\sigma, \mathbf{y}^\sigma)$$ for any permutation $\sigma\in\mathfrak{S}_n$. (Expressions like $\mathbf{x}^\sigma$ represent re-ordering the $x_i$ according to $\sigma$, so that $\mathbf{x}^\sigma = x_{\sigma(1)}, x_{\sigma(2)}, \ldots, x_{\sigma(n)}.$)
All these properties obviously hold for both $\operatorname{Var}$ and $C$ from inspecting the forms of expressions $(1)$ and $(2)$. The only one that might need any explanation is location invariance. However, a constant shift of values of the $x_i$ changes neither the residuals nor the differences:
$$x_i - \bar{x} = (x_i - a) - \overline{x - a}$$
and
$$x_j - x_i = (x_j - a) - (x_i - a).$$
Consequently, it is indeed obvious that the first version of $(1)$ and $(2)$ are location invariant.
Solution
Here, then, is the reasoning. Because $C$ is symmetric and multilinear, it is a quadratic form completely determined by coefficients $c_{ij} = c_{ji}$:
$$C(\mathbf{x}, \mathbf{y}) = \sum_{i, j=1}^n c_{ij}\, x_i y_j.$$
Because it is permutation-invariant, $c_{ij} = c_{i^\prime j^\prime}$ for any indices $i,j,i^\prime,j^\prime$ for which $i\ne j$ and $i^\prime \ne j^\prime$; also, $c_{ii} = c_{i^\prime i^\prime}$ for all indices $i$ and $i^\prime$. Thus, $C$ is determined by just two numbers, say $c_{11}$ and $c_{12}$. Finally, one of these numbers determines the other two by virtue of the location-invariance: that condition means
$$0 = C(\mathbf{0},\mathbf{0}) \overset{\text{location-invariance}}{=} C(\mathbf{1},\mathbf{0}) \overset{\text{symmetry}}{=} C(\mathbf{0},\mathbf{1}) \overset{\text{location-invariance}}{=} C(\mathbf{1},\mathbf{1})$$
(where "$\mathbf{0}$" and "$\mathbf{1}$" refer to constant $n$-vectors with these values). But
$$0=C(\mathbf{1},\mathbf{1}) = \sum_{i,j}^n c_{ij} = nc_{11} + (n^2-n)c_{12},$$ determining each of $c_{11}$ and $c_{12}$ in terms of the other.
This already proves the main point: $C$ must be proportional to $\operatorname{Cov}$, since each is determined by any single one of their coefficients. To find the constant of proportionality, inspect the two formulas $(1)$ and $(2)$, looking for all appearances of $x_1^2$: you can read off the associated value of $c_{11}$ from them. From the second version of $(1)$, the coefficient of $x_1^2$ clearly is $1/n - (1/n)^2$. From the first version of $(2)$, with $\mathbf{y} = \mathbf{x}$, the coefficient of $x_1^2$ clearly is $n-1$. (Geometrically, each point in the scatterplot of $(\mathbf{x},\mathbf{x})$ is paired with $n-1$ others, whence the square of its coordinate will appear $n-1$ times.) Therefore
$$c(n) = \frac{n-1}{1/n - (1/n)^2} = n^2,$$
QED. This was the only calculation required to demonstrate
$$\operatorname{Cov}(\mathbf{x}, \mathbf{y}) = \frac{1}{n^2}C(\mathbf{x}, \mathbf{y}) = \frac{1}{n^2}\sum_{i=1}^{n-1} \sum_{j=i+1}^n (x_j - x_i)(y_j - y_i).$$
