How is the relationship between two variables $X$ and $Y$ supposed to "explain" $R^2\text%$ of the variation of the data? Suppose we have a linear regression and we calculate $R^2 = 0.81$. What do we mean when we say "the relationship between two variables $X$ and $Y$ explains $81\text%$ of the variation of the data"? Especially this word "explains" confuses me, since the linear regression just reduces the variance in the observed data $Y$ by $81\text%$. It doesn't explain anything.

$R^2$ is a comparison of the variance of the error terms of two models. One model is naïve and always makes the same guess of $\bar y$ each time; this is how we get the "total sum of squares". The other model uses some features in an attempt to have less variance in the error term; this is how we get the "residual sum of squares".
In that regard, $R^2$ is a measure of the extent to which the error term variance is reduced. Some proportion of the original variance (from the naïve model) remains, but some proportion is explained by considering features.
Let's do an example. You want to know the height of a person you select from your family. Knowing nothing about the person you select, the best guess (in terms of minimizing square loss) is the mean height of people in your family; this is the naïve model. However, you know that adults tend to be taller than children. By considering age, you reduce the variability in the heights; this is the regression model of interest. Some of the reason the raw height values vary is because the subjects have different ages; age explains some of the variability in heights.
In math, and (shamelessly) taking some of a previous post of mine...
Notation
$y_i$ is observation $i$ of some response variable $Y$.
$\hat{y}_i$ is the value of $y_i$ predicted by the regression.
$\bar{y}$ is the average of all observations of the response variable.
$$ y_i-\bar{y} = (y_i - \hat{y_i} + \hat{y_i} - \bar{y}) = (y_i - \hat{y_i}) + (\hat{y_i} - \bar{y}) $$
$$( y_i-\bar{y})^2 = \Big[ (y_i - \hat{y_i}) + (\hat{y_i} - \bar{y}) \Big]^2  =  
(y_i - \hat{y_i})^2 + (\hat{y_i} - \bar{y})^2 + 2(y_i - \hat{y_i})(\hat{y_i} - \bar{y})
$$
$$SSTotal := \sum_i ( y_i-\bar{y})^2 = \sum_i(y_i - \hat{y_i})^2 + \sum_i(\hat{y_i} - \bar{y})^2 + 2\sum_i\Big[ (y_i - \hat{y_i})(\hat{y_i} - \bar{y}) \Big]$$
$$ :=SSRes + SSReg + Other $$
In OLS regression, $Other = 0$.
$$R^2= 1-\dfrac{SSRes}{SSTotal}\\
= \dfrac{SSTotal -SSRes-Other}{SSTotal} \\
= \dfrac{SSReg}{SSTotal} \\
=\dfrac{SSReg/n}{SSTotal/n} \\
=\dfrac{
\mathbb Var(\epsilon_{model})
}{
\mathbb Var(\epsilon_{naïve})
}$$
The claim about "proportion of variance explained" only applies in-sample (related to the $Other$ term not being zero), so if we consider the sample to be a population, then dividing by $n$ gives the exact variance.
