If $f$ is a strictly increasing real function such that $f(0)=0$ then using the Riemann series i need to calculate
$$\int_{[0,a]} f + \int_{[0,f(a)]} f^{-1}$$
From the fact that $f$ is increasing i can deduce that $f^{-1}$ is also increasing and therefore both are integrable but that is as far as i can get.
The Wiki link assumes that $f$ is continuous, and the other answers so far are not using Riemann sums as requested.
Let $0 = x_0 < x_1 < \cdots < x_n = a$, so $P_x = \{x_0,x_1,\ldots,x_n\}$ is a partition of $[0,a]$. Define $y_j = f(x_j)$ for each $0 \leq j \leq n$. Then $0 = f(0) = y_0 < y_1 < \cdots < y_n = f(a)$, and therefore $P_y = \{y_0,y_1,\ldots,y_n\}$ is a partition of $[0,f(a)]$.
Since $f$ is increasing, the upper and lower sums of $f$ with respect to $P_x$ are $$U(f,P_x) = \sum_{j=1}^{n} f(x_j)(x_j - x_{j-1})$$ and $$L(f,P_x) = \sum_{j=1}^{n} f(x_{j-1})(x_j - x_{j-1})$$ Since $f^{-1}$ is increasing, and $f^{-1}(y_j) = x_j$, and $f(x_j) = y_j$, the upper and lower sums of $f^{-1}$ with respect to $P_y$ are $$U(f^{-1},P_y) = \sum_{j=1}^{n} x_j(f(x_j) - f(x_{j-1}))$$ and $$L(f^{-1},P_y) = \sum_{j=1}^{n} x_{j-1}(f(x_j) - f(x_{j-1}))$$ Now, if we experiment with these sums a bit, we might notice that $$\begin{aligned} U(f,P_x) + L(f^{-1},P_y) &= \sum_{j=1}^{n}[f(x_j)(x_j - x_{j-1}) + x_{j-1}(f(x_j) - f(x_{j-1}))] \\ &= \sum_{j=1}^{n} [f(x_j)x_j - f(x_{j-1})x_{j-1}] \\ \end{aligned}$$ which telescopes to give us $f(x_n) x_n - f(0)x_0 = f(a)a - 0 = f(a)a$. This is rather nice because it gives us the exact answer even though we are only approximating with Riemann sums. However, note that this result depends on the special way we constructed $\{y_0,y_1,\ldots,y_n\}$ as $\{f(x_0),f(x_1),\ldots,f(x_n)\}$. Let's see if this is OK.
Note that both $f$ and $f^{-1}$ are Riemann integrable, because they are increasing functions. Let $\epsilon > 0$. Then there is a $\delta_x > 0$ such that $$0 \leq \int_0^a f(x)\ dx - U(f, P_x) < \epsilon$$ provided that all of the intervals in $P_x$ are less than $\delta_x$ in length. Similarly, there is a $\delta_y > 0$ such that $$-\epsilon < \int_0^{f(a)} f^{-1}(y)\ dy - L(f^{-1}, P_y) \leq 0$$ as long as all of the intervals in $P_y$ are less than $\delta_y$ in length. Therefore, as long as both interval length criteria are simultaneously met, we will be able to conclude that $$-\epsilon < \left(\int_0^a f(x)\ dx + \int_0^{f(a)} f^{-1}(y)\ dy\right) - (U(f, P_x) + L(f^{-1}, P_y)) < \epsilon$$ and as long as we maintain the condition $y_j = f(x_j)$, we will still have $U(f, P_x) + L(f^{-1}, P_y) = f(a)a$, and so $$\left|\int_0^a f(x)\ dx + \int_0^{f(a)} f^{-1}(y)\ dy - f(a)a\right| < \epsilon$$ So let's see how to satisfy the interval length criteria. Let us choose $\{x_0,x_1,\ldots,x_n\}$ freely such that $\max |x_j - x_{j-1}| < \delta_x$. This gives us a corresponding $\{y_0,y_1,\ldots,y_n\} = \{f(x_0),f(x_1),\ldots,f(x_n)\}$ which may or may not satisfy $\max |y_j - y_{j-1}| < \delta_y$. If not, then insert points $z_1,z_2,\ldots,z_m$ into $\{y_0,y_1,\ldots,y_n\}$ until the condition is satisfied, and then insert the corresponding $f^{-1}(z_1),f^{-1}(z_2),\ldots,f^{-1}(z_m)$ into $\{x_0,x_1,\ldots,x_n\}$, which only makes the latter partition finer, so it still satisfies its condition. Now both partitions are sufficiently fine, and they still satisfy $f(x_j) = y_j$ for every $j$, so we're done.