It's not much harder on the sphere than on the plane, once you recognize that
The points in question are the mutual intersections of three spheres: a sphere centered beneath location x1 (on the earth's surface) of a given radius, a sphere centered beneath location x2 (on the earth's surface) of a given radius, and the earth itself, which is a sphere centered at O = (0,0,0) of a given radius.
The intersection of each of the first two spheres with the earth's surface is a circle, which defines two planes. The mutual intersections of all three spheres therefore lies on the intersection of those two planes: a line.
Consequently, the problem is reduced to intersecting a line with a sphere, which is easy.
Here are the details. The inputs are points P1 = (lat1, lon1) and P2 = (lat2, lon2) on the earth's surface, considered as a sphere, and two corresponding radii r1 and r2.
Convert (lat, lon) to (x,y,z) geocentric coordinates. As usual, because we may choose units of measurement in which the earth has a unit radius,
x = cos(lon) cos(lat)
y = sin(lon) cos(lat)
z = sin(lat).
In the example, P1 = (-90.234036 Degree, 37.673442 Degree) has geocentric coordinates x1 = (-0.00323306, -0.7915, 0.61116) and P2 = (-90.953669 Degree, 36.109997 Degree) has geocentric coordinates x2 = (-0.0134464, -0.807775, 0.589337).
Convert the radii r1 and r2 (which are measured along the sphere) to angles along the sphere. By definition, one nautical mile (NM) is 1/60 degree of arc. Therefore, as angles,
r1 = 107.5 / 60 Degree = 0.0312705 radian
r2 = 145 / 60 Degree = 0.0421788 radian
The geodesic circle of radius r1 around x1 is the intersection of the earth's surface with a Euclidean sphere of radius sin(r1) centered at cos(r1)*x1.
The plane determined by the intersection of the sphere of radius sin(r1) around cos(r1)*x1 and the earth's surface is perpendicular to x1 and passes through the point cos(r1)*x1, whence its equation is x.x1 = cos(r1) (the "." represents the usual dot product); likewise for the other plane. There will be a unique point x0 on the intersection of those two planes that is a linear combination of x1 and x2. Writing x0 = a*x1 + b*x2 the two planar equations are
cos(r1) = x.x1 = (a*x1 + b*x2).x1 = a + b*(x2.x1)
cos(r2) = x.x2 = (a*x1 + b*x2).x2 = a*(x1.x2) + b
Using the fact that x2.x1 = x1.x2, which I shall write as q, the solution (if it exists) is given by
a = (cos(r1) - cos(r2)*q) / (1 - q^2),
b = (cos(r2) - cos(r1)*q) / (1 - q^2).
In the running example, I compute a = 0.973503 and b = 0.0260194.
Evidently we need q^2 != 1. This means that x1 and x2 can be neither the same point nor antipodal points.
Now all other points on the line of intersection of the two planes differ from x0 by some multiple of a vector n which is mutually perpendicular to both planes. The cross product
n = x1~Cross~x2
does the job provided n is nonzero: once again, this means that x1 and x2 are neither coincident nor diametrically opposite. (We need to take care to compute the cross product with high precision, because it involves subtractions with a lot of cancellation when x1 and x2 are close to each other.) In the example, n = (0.0272194, -0.00631254, -0.00803124).
Therefore, we seek up to two points of the form x0 + t*n which lie on the earth's surface: that is, their length equals 1. Equivalently, their squared length is 1:
1 = squared length = (x0 + t*n).(x0 + t*n) = x0.x0 + 2t*x0.n + t^2*n.n = x0.x0 + t^2*n.n
The term with x0.n disappears because x0 (being a linear combination of x1 and x2) is perpendicular to n. The two solutions easily are
t = sqrt((1 - x0.x0)/n.n)
and its negative. Once again high precision is called for, because when x1 and x2 are close, x0.x0 is very close to 1, leading to some loss of floating point precision. In the example, t = 1.07509 or t = -1.07509. The two points of intersection therefore equal
x0 + t*n = (0.0257661, -0.798332, 0.601666)
x0 - t*n = (-0.0327606, -0.784759, 0.618935)
Finally, we may convert these solutions back to (lat, lon) by converting geocentric (x,y,z) to geographic coordinates:
lon = ArcTan(x,y)
lat = ArcTan(Sqrt[x^2+y^2], z)
For the longitude, use the generalized arctangent returning values in the range -180 to 180 degrees (in computing applications, this function takes both x and y as arguments rather than just the ratio y/x; it is sometimes called "ATan2").
I obtain the two solutions (-88.151426, 36.989311) and (-92.390485, 38.238380), shown in the figure as yellow dots.
The axes display the geocentric (x,y,z) coordinates. The gray patch is the portion of the earth's surface from -95 to -87 degrees longitude, 33 to 40 degree latitude (marked off with a one degree graticule). The earth's surface has been made partly transparent to show all three spheres. The correctness of the computed solutions is evident by how the yellow points sit at the intersections of the spheres.