One geometric/probabilistic proof is by Eugenio Calabi.
Consider the double integral $I=\int_{0}^{1}\int_{0}^{1} \frac{1}{1-x^2y^2} \ dy \ dx.$ Expand the integrand into a geometric series as such $ \frac{1}{1-x^2y^2} =\sum_{n=0}^{\infty} (xy)^{2n},$ which is justified by the region of integration, that is $0 Exchange summation and integration (application of Tonelli's theorem), and integrate term by term to get $I=\sum_{n=0}^{\infty} \frac{1}{(2n+1)^2}.$ Now, we evaluate $I$ directly. Make the change of variables $x=\frac{\sin(u)}{\cos(v)}, y=\frac{\sin(v)}{\cos(u)}.$ This transformation has Jacobian Determinant $\frac{\partial (x,y) }{\partial(u,v)}=1-x^2y^2,$ which cancels with the integrand, and the region becomes the open triangle defined by the inequalities $u+v<\frac{\pi}{2}, u, v>0$ We know from geometry that this triangle has base $\frac{\pi}{2}$ and height $\frac{\pi}{2},$ hence the area is $\frac{\pi^2}{8}.$ Hence $I=\frac{\pi^2}{8}.$ Now to obtain $\zeta(2),$ we use the fact that $\zeta(2)=\sum_{n=1}^{\infty} \frac{1}{n^2} = \sum_{n=1}^{\infty} \frac{1}{(2n)^2}+\sum_{n=0}^{\infty} \frac{1}{(2n+1)^2},$ which implies by some algebraic manipulation that
$\zeta(2)=\frac{4}{3}I=\frac{\pi^2}{6}.$
We can extend this proof to the $\zeta(2k).$ Let
$I_{2k}= \int_{0}^{1}... \int_{0}^{1} \frac{1}{1-x_1^2 \dots x_{2k}^2} \ dx_{2k} \ ... \ dx_1.$
Converting this integrand into a geometric series and exchanging summation and integration gives the sum: $I_{2k}=\sum_{n=0}^{\infty} \frac{1}{(2n+1)^{2k}}.$ We can then recover $\zeta(2k)=\frac{2^{2k}}{2^{2k}-1} I_{2k},$ by observing $\zeta(2k)= \sum_{n=1}^{\infty} \frac{1}{(2n)^{2k}} + I_{2k},$ and using some algebraic manipulation.
On the other hand, we make the change of variables $x_i=\frac{\sin(u_i)}{\cos(u_{i+1})}, 1 \leq i \leq 2k,$ and $u_{2k+1}:=u_{1}$ here. This transformation turns out to have a Jacobian Determinant which cancels with the denominator of the integrand in $I_{2k}.$ The region of integration becomes the open polytope: $\Delta^{2k}= \left \lbrace (u_1, \dots, u_{2k}): u_{i}+u_{i+1} < \frac{\pi}{2}, u_i>0, 1 \leq i \leq 2k \right \rbrace$ in which $u_{2k+1}:=u_1.$ Computing the volume of this polytope is much more difficult. Rescaling with the change of variables $u_i=\frac{\pi}{2} v_i,$ and letting $V_1, \dots, V_{2k}$ being $2k$ independent, uniformly distributed random variables on $(0,1),$ we get
$I_{2k}=\text{Vol}(\Delta^{2k})=\left(\frac{\pi}{2}\right)^{2k} \text{Pr} \left( V_1+V_2<1, \dots, V_{2k-1}+V_{2k}<1, V_{2k}+V_{1}<1 \right),$ the probability that $V_1, \dots, V_{2k}$ have cyclically pairwise consecutive sums less than $1.$
In the literature, the aforementioned probability has been computed in several ways (see: https://www.maa.org/sites/default/files/pdf/news/Elkies.pdf https://arxiv.org/pdf/1003.3602.pdf https://pdfs.semanticscholar.org/35be/01e63c0bfd32b82c97d58ccc9c35471c3617.pdf)
Joshua Seaton also mentioned Zagier's and Kontsevich's reproduced Calabi proof. The general version of this is to evaluate
$J_{2k}= \int_{0}^{1}... \int_{0}^{1} \frac{1}{\sqrt{x_1 \dots x_{2k}}(1-x_1 \dots x_{2k})}\ dx_{2k} \ ... \ dx_1.$
A quick substitution shows $J_{2k}=2^{2k}I_{2k}.$ The generalized version of the change of variables is
$x_i=\frac{\xi_i^2(1+\xi_{i+1}^2)}{1+\xi_i^2}, \dots, x_{2k}=\frac{\xi_{2k}^2(1+\xi_{1}^2)}{1+\xi_{2k}^2},$ and upon computing its Jacobian Determinant, we get that
$J_{2k}=\int_{\mathbb{H}^{2k}} \frac{1}{\xi_1^2+1} \dots \frac{1}{\xi_{2k}^2+1} \ d \xi_{1} \dots d \xi_{2k},$ where $\mathbb{H}^{2k}$ is the hyperbolic defined by the inequalities: $\xi_{i}\xi_{i+1} <1, \xi_i>0,$ where $1 \leq i \leq 2k$ and $\xi_{2k+1}:=\xi_1.$ Letting $\Xi_1, \dots \Xi_{2k}$ be $2k$ independent, nonnegative Cauchy random variables, we get
$I_{2k}=\left(\frac{\pi}{2}\right)^{2k} \text{Pr} \left( \Xi_1\Xi_2<1, \dots, \Xi_{2k-1}\Xi_{2k}<1, \Xi_{2k}\Xi_{1}<1 \right),$ the probability that $\Xi_1, \dots, \Xi_{2k}$ have cyclically pairwise consecutive products less than $1.$
In my paper, I show that the probabilities lead to a very gargantuan formula:
$I_{2k}=\sum_{n=0}^{\infty} \frac{1}{(2n+1)^{2k}}= \left(\frac{\pi}{4} \right)^{2k} +\left(\frac{\pi}{4} \right)^{2k} \sum_{n=1}^{k} \sum_{\substack{(r_1, \dots, r_n) \in [2k]^n: \\ |r_p-r_q| \notin \lbrace 0,1,2k-1 \rbrace, \\ p,q \in [n]} } \prod_{i=1}^{n} \frac{1}{i+\sum_{j=1}^{i} \alpha_j},$ where $[m]:= \lbrace 1, \dots, m \rbrace$ and $\alpha_j=2- \delta(2k,2) - \sum_{m=1}^{j-1} \delta(|r_m-r_j|,2)+\delta(|r_m-r_j|,2k-2)$ and $\delta(a,b)$ is the Kronecker Delta Function. In particular, the inner sum is taken over all tuples $(r_1, \dots, r_n) \in [k]^ n$ having cyclically pairwise nonconsecutive entries.