If I use SnapPy to compute the Alexander polynomial of the link in the picture, I get $$t_1^2 - t_1 + 1$$ which is just the Alexander polynomail of the trefoil. But when I compute the Alexander polynomial of the knot exterior (the cusped 3-manifold),
L=Link([(1, 5, 2, 4), (3, 1, 4, 0), (7, 3, 0, 2), (5, 8, 6, 9), (9, 6, 8, 7)])
L.exterior().alexander_polynomial()
I get $$a^2 - ab + b^2$$ Is there a good explanation why the polynomials are different?

By a calculation through Fox calculus, the first elementary ideal of the presentation matrix for the infinite cyclic cover is $((1-t_1+t_1^2)(t_2-1),(1-t_1+t_1^2)(t_1-t_2))$, with $t_1$ corresponding to the trefoil and $t_2$ to the linked unknot. The GCD of this is the multivariable Alexander polynomial of a link, and this is $1-t_1+t_1^2$, matching the first output you say SnapPy produces. This presentation of the multivariable Alexander polynomial is special in that $t_1$ and $t_2$ are from the abelianizations of actual meridians.
The way SnapPy calculates the Alexander polynomial of a
Triangulation(the result of anexterior()of aLink) is to calculate the fundamental group and then going through the same Fox calculus computations I myself went through (seealexander_polynomial,alexander_polynomial_group, andalexander_polynomial_basicinpython/snap/nsagetools.py). However, when it computes the map to the abelianization of the fundamental group of the exterior, there is no guarantee that the algorithm will figure out the meridians of the link (and in general it cannot do so for general links, versus the case for knots due to Gordon and Luecke). The variables are only up to some isomorphism of the abelianization.The substitution $t_1=ab^{-1}$ and $t_2=b$ defines an isomorphism between $\mathbb{Z}[t_1^{\pm 1},t_2^{\pm 1}]$ and $\mathbb{Z}[a^{\pm 1},b^{\pm 1}]$. This carries the multivariable Alexander polynomial to $b^{-2}(a^2-ab+b^2)$. Since it is only defined up to multiplication by units in the group ring of the abelianization, this is equivalent to $a^2-ab+b^2$.