The question in this paper is stated in the bottom of page 28: "Suppose we are given , not just the maximum value for each year, but the largest ten (say) values. How might we use this data to obtain better estimates than could be made with just the annual maxima?".
With $r$ being the number of largest annual values, they are ordered (descending) and assumed to be iid. This is what enables eqn (2.1).
They don't just use the linear trend model introduced in equation (2.4) but also a quadratic model (eqn 3.1) and two different periodic models (11 and 19 years) (eqn 2.6).

Question: The authors want to compare which is the better estimate of the estimated quantity under different $r. $ How can I get
this conclusion by standard error?
For example, for MLE of $\alpha$, why is $r=5$ better than $r=1?$
Since we know that $$ \rm MSE=Bias^2+Variance, $$ do not we also need
the bias is smaller?
You can see in the comments here the reasons for using a biased estimator in order to get a lower MSE. Furthermore, with $r$ increasing from 1 to 5 or 10, we are bound to get a lower MSE: the $r$ observations are independent, hence the information function grows and we get a better Cramér–Rao lower bound.
Another question: for this non-stationary dataset, why did the author
set $\mu=\mu_n?$ This seems to be very similar to the residual in
linear regression model.
As you can see, they used different trend models. They do note just below eqn (3.3) the results obtained for the likelihood functions. You can see the comparison of the 4 models in fig. 2 where they compare the annual predicted maxima and medians. They further explain, under the summary section, that "In the Venice data, the main contribution to trend was a linear trend from 1931 onwards" and then explain that other trends had lesser effect.