You can employ the substitution $y=1- k^2 \sin^2x$ such that the integral can be written as
$$ K(k) = \int_{1-k^2}^1\!dy\,\frac{1}{2 \sqrt{y(1-y)(y-1 +k^2)}}.$$
Now, we can expand the integrand in $\epsilon = 1-k$ and obtain
$$ K(k) = \int_{1-k^2}^1\!dy\frac{1}{2 y \sqrt{1-y}} + O(1). \tag{1}$$
The estimate of the error term follows from the fact that expanding in $\eta =1-k^2$, we have
$$K(k) = \sum_{n=0}^\infty c_n \eta^n \int_{\eta}^1\!dy\, y^{-1-n} (1-y)^{-1/2}$$
with $c_n$ some constants. Now, we have for $n>0$ that
$$ \int_{\eta}^1\!dy\, y^{-1-n} (1-y)^{-1/2} = O(\eta^{-n})$$
such that only the $n=0$ term diverges for $\eta \to 0$ (which corresponds to $\epsilon \to 0$).
It thus remains to estimate the first term in (1) for $k \to 1$ which is not that difficult:
In fact due to the $1/y$ behavior close to $y=0$, the integral is logarithmically divergent and we have that
$$K(k) = \frac12 \left|\log (1-k^2)\right| + O(1) = \frac12 \left|\log \epsilon\right| + O(1).$$
Given the definite integral:
$$
I(a,\,b) := \int_0^{\frac{\pi}{2}} \frac{\cos\left(a\,\arcsin\left(b\,x\right)\right)}{\sqrt{1 - \left(b\,\sin x\right)^2}}\,\text{d}x
$$
through a substitution of the type $t = \frac{4}{\pi}\,x - 1$ we go back to this other integral:
$$
I(a,\,b) = \int_{-1}^1 \underbrace{\frac{\cos\left(a\,\arcsin\left(b\,\frac{\pi}{4}\,(t + 1)\right)\right)}{\sqrt{1 - \left(b\,\sin \left(\frac{\pi}{4}\,(t + 1)\right)\right)^2}}\,\frac{\pi}{4}}_{:= f(a,\,b,\,t)}\text{d}t
$$
to which it's possible to apply the Legendre-Gauss quadrature.
In particular, opting for the two-point formula:
$$
I(a,\,b) \approx k_1\,f(a,\,b,\,t_1) + k_2\,f(a,\,b,\,t_2)
$$
where is it:
$$
\begin{cases}
k_1\,t_1^0 + k_2\,t_2^0 = \frac{1 + (-1)^0}{1 + 0} \\
k_1\,t_1^1 + k_2\,t_2^1 = \frac{1 + (-1)^1}{1 + 1} \\
k_1\,t_1^2 + k_2\,t_2^2 = \frac{1 + (-1)^2}{1 + 2} \\
k_1\,t_1^3 + k_2\,t_2^3 = \frac{1 + (-1)^3}{1 + 3}
\end{cases}
\; \; \; \Leftrightarrow \; \; \;
\begin{cases}
k_{1,2} = 1 \\
t_{1,2} = \pm \frac{1}{\sqrt{3}}
\end{cases}
$$
it follows that:
$$
I(a,\,b) \approx f\left(a,\,b,\,-\frac{1}{\sqrt{3}}\right) + f\left(a,\,b,\,\frac{1}{\sqrt{3}}\right),
$$
approximation that involves at most an error equal to:
$$
\epsilon(a,\,b) = \frac{1}{135}\underset{-1 \le t \le 1}{\max} \left|\frac{\partial^4 f(a,\,b,\,t)}{\partial t^4}\right|.
$$
As an example:
$$
I\left(3,\,\frac{1}{10}\right) \approx f\left(3,\,\frac{1}{10},\,-\frac{1}{\sqrt{3}}\right) + f\left(3,\,\frac{1}{10},\,\frac{1}{\sqrt{3}}\right) = 1.51672
$$
with a maximum error equal to:
$$
\epsilon\left(3,\,\frac{1}{10}\right) = \frac{1}{135}\underset{-1 \le t \le 1}{\max} \left|\frac{\partial^4 f\left(3,\,\frac{1}{10},\,t\right)}{\partial t^4}\right| = 1.08341\cdot 10^{-4}\,.
$$
Best Answer
Since $\displaystyle K(k)=\int_0^{\pi/2}\frac{d t}{\sqrt{1-k^2\sin^2t}}$,
$$ \begin{align} &\frac{1}{\sqrt{(a+x)^2+z^2}}K\left(\frac{2\sqrt{ax}}{\sqrt{(a+x)^2+z^2}}\right) \\ =& \int_0^{\pi/2}\frac{dt}{\sqrt{(a+x)^2+z^2-4ax\sin^2t}} \\ =& \int_0^{\pi/2}\frac{dt}{\sqrt{a^2+z^2+x^2+2ax\cos2t}} \\=& \frac14\int_{-\pi}^{\pi}\frac{d\phi}{\sqrt{a^2+z^2+x^2+2ax\cos \phi}}\quad \phi=2t \end{align} $$
The integrand is the reciprocal of the distance between $(-x,\pi/2,0)$ and $\left(r=\sqrt{z^2+a^2},\theta=\arctan\dfrac az,\phi\right)$ under spherical coordinate, so one can expand it into series of spherical harmonics, assuming $x$ sufficiently small $$ \frac1{\sqrt{a^2+z^2+x^2+2ax\cos \phi}}=\sum_{l\ge0}\frac{(-x)^l}{r^{l+1}}\sum_{m=-l}^l(-1)^mP_l^{-m}(0)P_l^{m}(\cos\theta)e^{ i m\phi} $$ Perform the integral over $\phi$ and only the terms with $m=0$ remains, so $$ \int_{-\pi}^{\pi}\frac{d\phi}{\sqrt{a^2+z^2+x^2-2ax\cos \phi}}=2\pi\sum_{l\ge0}\frac{(-x)^l}{r^{l+1}}P_l(0)P_l(\cos\theta) $$ Now use Ramanujan's Master Theorem with $$ \varphi(l)=\frac\pi2\frac{1}{r^{l+1}}P_l(0)P_l(\cos\theta) $$
$$ \mathcal M\left(\frac{\pi}2\sum_{l\ge0}\frac{x^l}{r^{l+1}}P_l(0)P_l(\cos\theta)\right)(s)=\frac\pi{\sin\pi s}\varphi(-s)\\ =\Gamma(s)\Gamma(1-s)\cdot\frac\pi2\frac{1}{r^{1-s}}P_{-s}(0)P_{-s}(\cos\theta) $$
Finally, recall the special value of Legendre function $$ P_{\nu}(0)=\frac{\sqrt{\pi }}{\Gamma \left(\frac{1-\nu}{2}\right) \Gamma \left(\frac{\nu }{2}+1\right)} $$ and the duplication formula of Gamma function. After a little algebra, substitute$r=\sqrt{a^2+z^2},\cos\theta=\frac{z}{\sqrt{a^2+z^2}}$ and the result follows.