The following interesting derivation is from Aigner's "A Course in Enumeration" (Springer, 2007).
Remember that if we have the following exponential generating functions:
$\begin{align}
\widehat{A}(z)
&= \sum_{n \ge 0} a_n \frac{z^n}{n!} \\
\widehat{B}(z)
&= \sum_{n \ge 0} b_n \frac{z^n}{n!}
\end{align}$
then also:
$\begin{align}
\widehat{A}(z) \cdot \widehat{B}(z)
&= \sum_{n \ge 0}
\left(
\sum_{0 \le k \le n} \frac{a_k}{k!} \frac{b_{n - k}}{(n - k)!}
\right) z^n \\
&= \sum_{n \ge 0}
\left(
\sum_{0 \le k \le n} \frac{n!}{k! (n - k)!} a_k b_{n - k}
\right) \frac{z^n}{n!} \\
&= \sum_{n \ge 0}
\left(
\sum_{0 \le k \le n} \binom{n}{k} a_k b_{n - k}
\right) \frac{z^n}{n!}
\end{align}$
Let's define:
$\begin{align}
S_m(n)
= \sum_{1 \le k \le n - 1} k^m
\end{align}$
We can define the exponential generating function:
$\begin{align}
\widehat{S}_n(z)
&= \sum_{m \ge 0} S_m(n) \frac{z^m}{m!} \\
&= \sum_{1 \le k \le n - 1} \sum_{m \ge 0} \frac{k^m z^m}{m!} \\
&= \sum_{1 \le k \le n - 1} \mathrm{e}^{k z} \\
&= \frac{\mathrm{e}^{n z} - 1}{\mathrm{e}^z - 1} \\
\end{align}$
This is almost the exponential generating function of the powers of $n$:
$\begin{align}
\widehat{P}_n(z)
&= \sum_{m \ge 0} n^m \frac{z^m}{m!} \\
&= \mathrm{e}^{n z}
\end{align}$
We can write:
$\begin{align}
(\widehat{P}_n(z) - 1) \widehat{B}(z)
= z \widehat{S}_n(z) \tag{1}
\end{align}$
where we have the exponential generating function of the Bernoulli numbers:
$\begin{align}
\widehat{B}(z)
&= \frac{z}{\mathrm{e}^z - 1} \\
&= \sum_{n \ge 0} B_n \frac{z^n}{n!}
\end{align}$
Comparing coefficients of $z^{m + 1}$ in (1):
$\begin{align}
\sum_{m \ge 1} z^m
\sum_{0 \le k \le m}
\binom{m}{k} \frac{(n z)^{m - k}}{(m - k)!} B_k
= \sum_{m \ge 0} S_m(n) \frac{z^{m + 1}}{m!}
\end{align}$
we get after simpĺifying:
$\begin{align}
S_m(n)
= \frac{1}{m + 1} \,
\sum_{0 \le k \le m} \binom{m + 1}{k} B_k n^{m + 1 - k}
\end{align}$
Note: The formula given is often associated with Faulhaber, but Faulhaber's formulas where quite different (and computationally more efficient). This formula is due to Bernoulli.
Since you asked for an intuitive explanation consider a simple case of $1^2+2^2+3^2+4^2$ using a set of children's blocks to build a pyramid-like structure.
First you arrange $16$ blocks in a $4\times4$ square. Next you place on top of this arrangement, $9$ blocks in a $3\times3$ square aligning the blocks of the upper left corner of each square, one above the other. On top of this construction place $4$ blocks in a $2\times2$ square, similarly aligned. Finally crown the upper left corner with a single block for the $1\times1$ square.
The following table represents the number of block in each column of the arrangement viewed from above:
$\begin{array}{cccc}
4&3&2&1\\
3&3&2&1\\
2&2&2&1\\
1&1&1&1
\end{array}$
We can find the total number of blocks in the arrangement by adding the number of columns containing a single block to twice the number of columns containing two blocks then adding three times the number of columns containing three blocks and finally adding four times the one column containing four blocks.
\begin{eqnarray}
\sum_{i=1}^4i^2=1\cdot(4+3)+2\cdot(3+2)+3\cdot(2+1)+4\cdot1
\end{eqnarray}
If we do the same thing with a pyramid of blocks having $n\times n$ blocks on its base then the summation would look like the following:
\begin{eqnarray}
\sum_{i=1}^{n}i^2&=&1\cdot(n+n-1)+2\cdot(n-1+n-2)+3\cdot(n-2+n-3)\\
&+&\cdots+(n-1)\cdot(2+1)+n\cdot1\\
&=&1\cdot(2n-1)+2\cdot(2n-3)+3\cdot(2n-5)+\cdots+(n-1)\cdot3+n\cdot1\\
&=&\sum_{i=1}^{n}i(2n-2i+1)\\
&=&2n\sum_{i=1}^{n}i-2\sum_{i=1}^{n}i^2+\sum_{i=1}^{n}i\\
\sum_{i=1}^{n}i^2&=&(2n+1)\sum_{i=1}^{n}i-2\sum_{i=1}^{n}i^2\\
3\sum_{i=1}^{n}i^2&=&(2n+1)\sum_{i=1}^{n}i\\
&=&\dfrac{n(n+1)(2n+1)}{2}\\
\sum_{i=1}^{n}i^2&=&\dfrac{n(n+1)(2n+1)}{6}
\end{eqnarray}
Best Answer
The binary expansion of $\sum_{k=0}^n2^k$ is a string of $n+1$ 1's: $$\underbrace{111\dots111}_{n+1}$$ If I add a 1 to this number, what do I get? $$1\underbrace{000\dots000}_{n+1}$$ 1 followed by $n+1$ 0's, hence $2^{n+1}$. Therefore $$\sum_{k=0}^n2^k=2^{n+1}-1$$