Multilateration of a light source in a medium whose index of refraction varies with position

Background

I have a set of $N$ receivers whose locations in $3D$ space are well-known. These receivers are immersed in a medium whose index of refraction (and thereby the velocity of light propagation) varies as a function of depth (namely, as $n(z) = A - B e^{Cz}$ for some positively-valued constants $A,B$ and $C$). A light source whose location is unknown emits a signal at some unknown time. Each of the $N$ receivers records the time of signal arrival. I wish to use this data to reconstruct the $3D$ coordinates of the unknown light source.

What I've Tried

I have experience solving this problem in the simpler case wherein the velocity of the signal is constant throughout the medium. I've written or employed simple least squares minimization algorithms, the Bancroft algebraic solution, the Chan algebraic solution, the Taylor algorithm, and the Taylor-Chan algorithm.

Regrettably, these methods work quite poorly in this more complex case, and try as I may I cannot seem to find an algebraic nor iterative numerical adaptation or solution which is effective. The range of possible coordinates is sufficiently large that "stupid simple" solutions such as simply iterating very finely over some domain are far too slow and computationally expensive for my use case, as are e.g. simple "approximation search" algorithms and the like.

Some Observations & Parameters

1. I know the receivers and source to be located underground (i.e. $z < 0$), which may be useful in some algebraic solution wherein there exist multiple solutions.
2. The index of refraction is a function only of depth, and I (well, mostly some other folks) have solved the appropriate Euler-Lagrange equations to obtain an analytical expression for the path taken by a light ray between two points.
3. I'd prefer to employ as few receivers as possible, however I can use as many as $20$.
4. I'd prefer an algorithm wherein I don't have to worry too much about adjusting parameters for use with e.g. specific receiver configurations.

What else might I try? What solutions/research exists on this topic, if any? I've tried perusing the literature, however much of it is couched in terms of terminology that is slightly unfamiliar to me as someone with a background in physics as opposed to something like applied mathematics or signal processing. I am unsure of what to search for, and expect that this is surely something that has already been studied.

There's a similar (but slightly different) question here, whose answers unfortunately aren't terribly useful in my case.

Edit: A Possible Lead

Given points $\langle x_0, y_0, z_0 \rangle$ and $\langle x_1, y_1, z_1 \rangle$, we can compute the transit time for a ray propagating between these points as:

$$ A\sqrt{\frac{C^{2}+K^{2}}{C^{2}K^{2}}}\left( \ln\left(\left|\frac{\sqrt{\frac{A+Be^{Cz_{1}}-\beta}{A+Be^{Cz_{1}}+\beta}}-\alpha}{\sqrt{\frac{A+Be^{Cz_{1}}-\beta}{A+Be^{Cz_{1}}+\beta}}+\alpha}\right|\right) - \ln\left(\left|\frac{\sqrt{\frac{A+Be^{Cz_{0}}-\beta}{A+Be^{Cz_{0}}+\beta}}-\alpha}{\sqrt{\frac{A+Be^{Cz_{0}}-\beta}{A+Be^{Cz_{0}}+\beta}}+\alpha}\right|\right)\right) $$

—where:

$$ K = C\frac{\sqrt{A^2-\beta^2}}{\beta} \qquad\textbf{ and }\qquad \alpha=\sqrt{\frac{A-\beta}{A+\beta}} $$

—and where $\beta$ denotes the angle made between the light ray and the horizontal axis at the initial point $\langle x_0, y_0, z_0 \rangle$.