You're surely right that there cannot be a "closed form" for such a series;
but it can still be approximated to any desired precision.
The defining sum
$$
x = x(a)
= \sum_{i=0}^\infty \sum_{j=0}^\infty \exp\Bigl(-a \sqrt{i^2+j^2}\Bigr)
$$
requires only $O(N^2/a)$ terms to approximate to within $\exp(-N)$,
so it converges quickly enough unless $a$ is small;
but in that case one can use a Laurent series
$x(a) = \frac\pi{2a^2} + \sum_{l=0}^\infty c_l a^l$
that converges for $a < 2\pi$. The "small" and "large" ranges overlap,
so one can check the analysis by computing $x(a)$ both ways
for some values of $a$ in say $[1/2, 2]$ and noting that
they agree to within the computed precision.
Instead of $x(a)$, which is a sum over integer $(i,j)$
in the first quadrant, it is more natural to sum over all integer pairs:
define
$$
F(a) = \sum_{i=-\infty}^\infty \sum_{j=-\infty}^\infty
\exp \Bigl(-a \sqrt{i^2+j^2}\Bigr).
$$
It is easy to go between $x(a)$ and $F(a)$,
because each term in the $x(a)$ sum appears $4$ times in the $F(a)$ sum
(for $(i,j)$, $(i,-j)$, $(-i,j)$, and $(-i,-j)$),
except when $i=0$ or $j=0$ --- and those terms do give geometric series,
so we soon find
$$
F(a) = 4 x(a) - \frac4{1 - e^{-a}} + 1.
$$
We shall prove that for $0 < a < 2\pi$ we have
$$
F(a) = \frac{2\pi}{a^2} + \sum_{k=1}^\infty c_{2k-1} a^{2k-1},
$$
where the coefficients $c_1,c_3,c_5,\ldots$ are given by
$$
c_{2k-1} = 4 \left({-3/2 \atop k-1}\right) (2\pi)^{-2k}
\zeta\bigl(k+\frac12\bigr) L\bigl(k+\frac12, \chi_4\bigr).
$$
Here $\zeta$ is the Riemann zeta function and
$L(\cdot, \chi_4)$ is the Dirichlet $L$-function
$$
L(s, \chi_4) = \sum_{n=0}^\infty \frac{(-1)^n}{(2n+1)^s}
= 1 - \frac1{3^s} + \frac1{5^s} - \frac1{7^s} + - \cdots.
$$
These factors of $c_{2k-1}$ are not elementary but can be
efficiently approximated to any desired accuracy,
so the same is true of the Laurent series. Here is a gp implementation:
L4(s) = sumalt(n=0,(-1)^n/(2*n+1)^s)
{
F(a) = 2*Pi/a^2 + 8*Pi * suminf(k=1,
binomial(-3/2,k-1) * (2*Pi)^(-2*k-1) * zeta(k+1/2) * L4(k+1/2) * a^(2*k-1))
}
To obtain the series expansion of $F(a)$, we first apply
the 2-dimensional
Poisson
summation formula
$$
\sum_{i=-\infty}^\infty \sum_{j=-\infty}^\infty f(i,j)
= \sum_{m=-\infty}^\infty \sum_{n=-\infty}^\infty \hat f(m,n),
$$
where $\,f: {\bf R}^2 \to \bf C$ is "any" function, and
$\,\hat f$ is its Fourier transform
$$
\hat f(y_1,y_2) = \mathop{\iint}_{(x_1,\,x_2) \in {\bf R}^2}
e^{-2\pi i (x_1 y_1 + x_2 y_2)} f(x_1,x_2) \, dx_1 \, dx_2.
$$
For $f(x_1,x_2) = \exp\bigl( -a \sqrt{x_1^2+x_2^2} \bigr)$
the Fourier transform is elementary, though it takes a bit of work
(write the double integral in polar coordinates,
and integrate with respect to $r$ first). We find
$$
F(a) = 2 \pi a \sum_{m=-\infty}^\infty \sum_{n=-\infty}^\infty
\bigl(a^2 + 4\pi^2(m^2+n^2)\bigr)^{-3/2}.
$$
This sum converges, but too slowly to use as it stands.
(The singularity of $\,f$ at $(0,0)$ makes $\,\hat f$ decay slowly
at infinity, though the decay is still fast enough to justify
the use of Poisson summation.) So for each $(m,n) \neq (0,0)$
we write $(a^2 + 4\pi^2(m^2+n^2))^{-3/2}$ as a binomial expansion
$$
(2\pi)^{-3} (m^2+n^2)^{-3/2} \sum_{k=1}^\infty
\left({-3/2 \atop k-1}\right) \bigl(a^2 / (4\pi^2 (m^2+n^2))\bigr)^{k-1},
$$
which converges if $|a| < 2 \pi$, and then combine terms with
$a^2$ appearing to the same power. The formula for the coefficients
$c_{2k-1}$ then follows using the known factorization
$$
\mathop{\sum\!\sum}_{\scriptstyle(m,n) \in {\bf Z}^2 \atop \scriptstyle(m,n) \neq (0,0)}
(m^2+n^2)^{-s} = 4 \zeta(s) \, L(s,\chi_4)
$$
which encodes unique factorization in ${\bf Z}[i]$.
