Let $0\leq i\leq k$ and $\alpha,\beta>0$ such that the image of $f\big|_{(-\alpha,\alpha)\times(t_{i}-\beta,t_{i}+\beta)}$ is contained in $U\subseteq M$ where $(U,\rho)$ is a coordinate neighbourhood of $f(0,t_{i})$. Then there exist continuous $\lambda_{1},\ldots,\lambda_{n}:(-\alpha,\alpha)\times(t_{i}-\beta,t_{i}+\beta)\to\mathbb{R}$ such that $$\frac{\partial f}{\partial s}(s,t)=\sum_{j=1}^{n}\lambda_{j}(s,t)\frac{\partial}{\partial x_{j}}\bigg|_{f(s,t)}$$ It follows from the definition of the variation that $\lambda_{j}\big|_{(-\alpha,\alpha)\times[t_{i},t_{i}+\beta)}$ and $\lambda_{j}\big|_{(-\alpha,\alpha)\times(t_{i}-\beta,t_{i}]}$ are differentiable. Furthermore it is a defining property of the covariant derivation that if well defined, we obtain:
$$\frac{D}{\partial s}\bigg|_{(s,t)}\frac{\partial f}{\partial s}=\sum_{j=1}^{n}\frac{\partial \lambda_{j}}{\partial s}(s,t)\frac{\partial}{\partial x_{j}}\bigg|_{f(s,t)}+\sum_{j=1}^{n}\lambda_{j}(s,t)\frac{D}{\partial s}\bigg|_{(s,t)}\frac{\partial}{\partial x_{j}}\\
=\sum_{j=1}^{n}\frac{\partial \lambda_{j}}{\partial s}(s,t)\frac{\partial}{\partial x_{j}}\bigg|_{f(s,t)}+\sum_{j,m=1}^{n}\lambda_{j}(s,t)\lambda_{m}(s,t)\left(\nabla_{\frac{\partial}{\partial x_{m}}}\frac{\partial}{\partial x_{j}}\right)_{f(s,t)}$$
So all there is to check is that $\frac{\partial \lambda_{j}}{\partial s}(s,t_{i}^{+})=\frac{\partial \lambda_{j}}{\partial s}(s,t_{i}^{-})$. But this follows from the definition of $\frac{\partial \lambda_{j}}{\partial s}(s,t_{i})$. More formally, we can define $\lambda_{j}^{1}:=\lambda_{j}\big|_{(-\alpha,\alpha)\times(t_{i}-\delta,t_{i}]}$ and $\lambda_{j}^{2}:=\lambda_{j}\big|_{(-\alpha,\alpha)\times[t_{i},t_{i}+\delta)}$. By definition we have:
$$\frac{\partial \lambda_{j}}{\partial s}(s,t_{i}^{+})=\lim_{t\downarrow t_{i}}\frac{\partial\lambda_{j}}{\partial s}(s,t)=\lim_{t\downarrow t_{i}}\frac{\partial\lambda_{j}^{2}}{\partial s}(s,t)$$
and:
$$\frac{\partial \lambda}{\partial s}(s,t_{i}^{-})=\lim_{t\uparrow t_{i}}\frac{\partial\lambda_{j}}{\partial s}(s,t)=\lim_{t\uparrow t_{i}}\frac{\partial\lambda_{j}^{1}}{\partial s}(s,t)$$
As the restrictions are continuously differentiable by assumption, we obtain $\frac{\partial \lambda_{j}}{\partial s}(s,t_{i}^{+})=\frac{\partial \lambda_{j}^{2}}{\partial s}(s,t_{i})$ and $\frac{\partial \lambda_{j}}{\partial s}(s,t_{i}^{-})=\frac{\partial \lambda_{j}^{1}}{\partial s}(s,t_{i})$. Explicit calculation now yields:
$$\frac{\partial\lambda_{j}}{\partial s}(s,t_{j}^{+})=\lim_{h\to 0}\frac{\lambda_{j}^{2}(s+h,t_{i})-\lambda_{j}^{2}(s,t_{i})}{h}=\lim_{h\to 0}\frac{\lambda_{j}(s+h,t_{i})-\lambda_{j}(s,t_{i})}{h}\\ \frac{\partial\lambda_{j}}{\partial s}(s,t_{j}^{-})=\lim_{h\to 0}\frac{\lambda_{j}^{1}(s+h,t_{i})-\lambda_{j}^{1}(s,t_{i})}{h}=\lim_{h\to 0}\frac{\lambda_{j}(s+h,t_{i})-\lambda_{j}(s,t_{i})}{h}$$
So indeed this is independent of whether we take the derivative from below $t_{i}$ or from above $t_{i}$.
As a synopsis: as expected this is nothing but the fact that at $t_{i}$ the definition of a variation leads to the smooth curve $s\mapsto f(s,t_{i})$. Hence all your favourite operations can be applied to this curve and as $t\mapsto f(s,t)$ is continuous, also the derivatives are continuous in the parameter.
P.s.:@Jason: I am sure that this just boils down to the solutions of a differential equation being continuous with respect to the initial conditions. The formal discussion just covers a special case and indeed the proof is nothing but the fact that continuity means continuity from both sides. Thank you for your patience anyway.
This is a rough sketch of the proof, I still have to check for correctness.
We have $\exp_{y_s} (tX(s))$ as a coordinate-field, not a vector field. However, when we do operations like $\frac{D}{ds}$, they correspond to covariant derivatives along the directions of the field given by
$\frac{D}{ds}\exp_{y_s} (tX(s))$.
$$\frac{D^2}{dt} J(t) = \frac{D}{dt} \frac{D}{dt} (\frac{D}{ds} \exp_{y_s} (tX(s))) = \\= \frac{D}{dt} \frac{D}{ds} \frac{D}{dt}\exp_{y_s} (tX(s)) = \frac{D}{dt} \frac{D}{ds} \dot \tau(t) = \\=\frac{D}{ds} \frac{D}{dt} \dot \tau(t) \; - R(ds,dt)\dot \tau(t) =\\= -R(J(t), \dot \tau(t))\dot \tau(t) $$
We can commute the two innermost covariant derivatives because we start off with a coordinate field. (it works roughly for the same reason that $s_{;ab} = s_{;ba}$ for a scalar field). For commuting the two outermost covariant derivatives we apply the Ricci identity. Also $\frac{D}{dt} \dot \tau(t)$ vanishes because they are geodesics.
The Jacobi equation is a second order differential equation whose solutions are completely determined by the values of $J(p)$ and $\dot J(p)$ for a given parameter $p$ along the geodesic. At $\tau(0) = x_0$ we have $J(0) = V$ and
$$\dot J(0) = \frac{D}{dt} \frac{D}{ds} \exp_{y_s} (tX(s))) = \\= \frac{D}{ds} \frac{D}{dt} \exp_{y_s} (tX(s))) =
\frac{D}{ds} \dot \tau(t)$$
at values $s= 0,t=0$. Given that the initial geodesic family is formed by parallel-transporting $\dot \tau(0)$ along the geodesic spanned by $ds$, we get that $\frac{D}{ds} \dot \tau(t) = 0$ at $s=0,t=0$ . Thus our initial conditions are $J(0) = V$, $\dot J(0) = 0$, which, along with the Jacobi equation, uniquely determine the field along the geodesic.
Best Answer
Because $\dfrac{\partial\phi}{\partial s}(0,a) = \dfrac{\partial\phi}{\partial s}(0,b) = 0$! As you yourself wrote, the endpoints stay fixed during the variation. Indeed, $\dfrac{\partial\phi}{\partial s}(s,a) = \dfrac{\partial\phi}{\partial s}(s,b) = 0$ for all $s$.