Apologies in advance for the long question.
My understanding is that in GR, massive observers move along timelike curves $x^\mu(\lambda)$, and if an observer moves from point $x^\mu(\lambda_a)$ to $x^\mu(\lambda_b)$, then his clock will measure that an amount of time $t_{ba}$ given by the curve's arc length; $$ t_{ba} = \int_{\lambda_a}^{\lambda_b}d\lambda \sqrt{-g_{\mu\nu}(x(\lambda))\dot x^\mu(\lambda)\dot x^\nu(\lambda)} $$ will have elapsed where $g_{\mu\nu}$ is a metric on spacetime with signature $(-,+,+,+)$.
Why is this so?
Here is how I would attempt to justify this fact in special relativity with $g_{\mu\nu} = \eta_{\mu\nu}$. Consider an inertial observer $O$ in $\mathbb R^{3,1}$, and suppose that this observer sees a clock, which I'll call observer $O'$ moving around on a curve $x^\mu(\lambda)$. If $O'$ were also an inertial observer, then given any event with coordinates $x^\mu$ as measured by $O$, observer $O'$ would measure the coordinates of the event to be $x'^\mu = \Lambda^\mu_{\phantom\mu\nu} x^\nu + x_0^\mu$ for some Lorentz transformation $\Lambda$. If $O'$ is not inertial, then this is no longer true, and there is some more complicated family of transformations, say $T_\lambda$ between events as seen by both observers.
I would argue, however, that if we were to partition the interval $[\lambda_a, \lambda_b]$ into a large number $N$ of intervals $I_1=[\lambda_a, \lambda_i], \dots, I_N=[\lambda_{N-1}, \lambda_b]$ with $\lambda_n = \lambda_a+n\epsilon_N$ and $\epsilon_N=(\lambda_b-\lambda_a)/N$, then on each interval $I_n$, $O'$ is approximately an inertial observer in the sense that $$ T_{\lambda_n} = P_n + \mathcal O(\epsilon_N), \qquad (\star) $$ for some Poincare transformation $P_n$. Then we would note that since $O'$ is stationary in his own reference frame, he measures his worldline to have the property $\dot x'^\mu(\lambda) = (\dot t(\lambda), \mathbf 0)$ so that $$ I_{ba}=\int_{\lambda_a}^{\lambda_b}d\lambda \,\sqrt{-\eta_{\mu\nu}\dot x'^\mu\dot x'^\mu} = \int_{\lambda_a}^{\lambda_b} d\lambda \, \sqrt{\dot t^2} = t(\lambda_b) - t(\lambda_a) = t_{ba} $$ On the other hand the integral on the left can be written as a Riemann sum using the partition above, and we can invoke ($\star$) above to get \begin{align} I_{ba} &= \lim_{N\to\infty}\left[\sum_{n=1}^N \epsilon_N\sqrt{-\eta_{\mu\nu}\dot x'^\mu(\lambda_n)\dot x'^\nu(\lambda_n)}\right] \notag\\ &= \lim_{N\to\infty}\left[\sum_{n=1}^N \epsilon_N\sqrt{-\eta_{\mu\nu}\dot x^\mu(\lambda_n)\dot x^\nu(\lambda_n)} + \mathcal{O}(\epsilon_N^2)\right] \notag\\ &= \int_{\lambda_a}^{\lambda_b}d\lambda \,\sqrt{-\eta_{\mu\nu}\dot x^\mu\dot x^\mu} \end{align} Combining these two computations gives the desired result.
How do others feel about this argument?
I'm not completely comfortable with it because of the assumption $(\star)$ I made on $T_\lambda$.
I imagine that in GR a similar argument could be made by invoking local flatness of the metric.