Jelly, 25 bytes

I;°ḊÆS,ÆẠƲP€SÆAHÆSḤ,Ʋ×⁽ßx

A dyadic Link that accepts the longitudes on the left and the latitudes, in the same order, on the right and yields the two distances, [direct, surface] in kilometres.

Try it online!

How?

Labelling the two longitudes as \$\lambda_i\$ and the two latitudes as \$\phi_i\$, the angle between the lines from each of the two points and the centre of the sphere, the central angle, is^[1]:

$$\Delta\sigma = \arccos{(\cos{|\lambda_1-\lambda_2|}\cos{\phi_1}\cos{\phi_2}+\sin{\phi_1}\sin{\phi_2})}$$

Note: Since \$\cos\$ is symmetric around the y-axis, we can use \$\lambda_2-\lambda_1\$ in place of \$|\lambda_1-\lambda_2|\$.

The shortest distance between two points in Euclidean space is a straight line, so we have an isosceles triangle with the equal side lengths being the sphere's radius, \$r=6371\$ kilometres, so the direct path in space is (using \$\text{opposite}=\text{hypotenuse}\sin\alpha\$):

$$d = 2r \sin \frac{\Delta\sigma}{2}$$

The shortest distance between two points on a sphere is the minor arc of the great circle on which they lie. This is the arc projected by the angle \$\Delta\sigma\$ and thus has length:

$$s = r \Delta\sigma$$

The code takes \$(\lambda_1, \lambda_2)\$ on the left and \$(\phi_1, \phi_2)\$ on the right and produces \$(d, s)\$:

I;°ḊÆS,ÆẠƲP€SÆAHÆSḤ,Ʋ×⁽ßx - Link: Longitudes; Latitudes
I                         - forward difference {Longitudes} -> [l2 - l1]
 ;                        - concatenate {Latitudes} -> [l2 - l1, p1, p2]
  °                       - convert from degrees to radians
         Ʋ                - last four links as a monad f(X=that):
   Ḋ                      -   dequeue {X} -> [p1, p2]
    ÆS                    -   sine -> [sin(p1), sin(p2)]
       ÆẠ                 -   cosine {X} -> [cos(l2 - l1), cos(p1), cos(p2)]
      ,                   -   pair -> [[sin(p1), sin(p2)],
                                       [cos(l2 - l1), cos(p1), cos(p2)]]
          P€              - product of each
            S             - sum
             ÆA           - arccosine -> CentralAngle
                    Ʋ     - last four links as a monad - f(CentralAngle):
               H          -   halve
                ÆS        -   sine
                  Ḥ       -   double
                   ,      -   pair -> [2 × sin(CentralAngle / 2),
                                       CentralAngle]
                      ⁽ßx - 6371
                     ×    - multiply -> [6371 × 2 × sin(CentralAngle / 2),
                                         6371 × CentralAngle]
                                      = [direct, surface]

Ruby, 105 bytes

->z,p,q{a,b,c,d,e=(z<<p-q).flat_map{|j|(1i**j/=90r).rect}
p Math.acos(m=a*c*e+b*d)*r=6371,(2-m*2)**0.5*r}

4 bytes saved thanks to Dingus.

Try it online!

Ruby, 109 bytes

A function taking 3 arguments: an array of latitudes, followed by the longitudes as separate arguments.

->z,p,q{a,b,c,d,e=(z<<p-q).map{|j|(1i**j/=90.0).rect}.flatten
p Math.acos(m=a*c*e+b*d)*r=6371,(2-m*2)**0.5*r}

Try it online!

It uses essentially the same formula as Jonathan Allen's answer, which can be found at Haversine formula . Ruby syntax for trig functions requires either include Math at the beginning of the program or Math. before every instance of a trig function. Therefore extensive use is made of complex numbers, which are available by default. The formula works as follows.

The cosine of the angle between the two points is calculated as a dot product (for a pair of vectors of magnitude 1, the dot product is equal to cosine of the angle between them):

u1*u2 + v1*v2 + z1*z2

In order to do this the latitude information for each point must first be converted from an angle into vectors r1+z1i and r2+z2i.

The longitudes may be converted to u+vi coordinates r1*(x1+y1i) and r2*(x2+y2i) but by rotating the coordinate system (subtracting longitude 1 from both longitudes 1 and 2) longitude 1 becomes effectively zero, hence there is no y contribution and the x contribution becomes 1. The dot product is therefore

cos theta = (r1*1)*(r2*x2) + 0 + z1*z2 or equivalently a*1*c*e + b*d

To avoid using the longwinded Math.sin function, the sine of half the angle can be found using the identity sin(theta/2) = (1/2-(cos theta)/2)**0.5. Because we want to know double this value, the code considers coefficients 4 times larger (2-m*2) to give the required value directly.

Commented code

->z,p,q{                        #Take input as an array of latitudes, plus longitudes as separate arguments
a,b,c,d,e=                      #Set up 5 variables for the complex numbers (the 6th will be discarded)
(z<<p-q).map{|j|                #Calculate the difference between longitudes and append to array. Iterate through 3 elements
(1i**j/=90.0).rect              #For each element convert angle into rectangular coordinates by raising i to a power
}.flatten                       #Flatten the array of 3 pairs and assign to a,b,c,d,e
p Math.acos(m=a*c*e+b*d)*r=6371,#Calculate the cosine of the angle and the surface distance
(2-m*2)**0.5*r}                 #Convert the cosine into the straight line distance.

Charcoal, 66 bytes

⊞θ⁻⊟θ⊟θ≔ΣＥ⪪.imag.real⁵ΠＥ…θ⁺²κＵＶ⁺ＩＸＩ1j∕λ⁹⁰ιθＩ×⁶³⁷¹⟦▷math.acosθ₂⁻²⊗θ

Attempt This Online! Link is to verbose version of code. Takes input as an array of two latitudes and two longitudes. Explanation: Port of @LevelRiverSt's Ruby answer.

⊞θ⁻⊟θ⊟θ

Subtract one longitude from the other.

≔ΣＥ⪪.imag.real⁵ΠＥ…θ⁺²κＵＶ⁺ＩＸＩ1j∕λ⁹⁰ιθ

Obtain the sine and cosine of each angle by dividing by 90, raising i to the resulting power, and extracting the real or imaginary part as appropriate. Multiply all three cosines and the sines of the latitudes together, then take the sum. This results in the cosine of the angle between the two points.

Ｉ×⁶³⁷¹⟦▷math.acosθ₂⁻²⊗θ

Take the arccosine to obtain the angle and also subtract the the doubled cosine from 2 and take the square root of the result. Multiply both values by the Earth's radius to obtain the desired output.

Would be 56 bytes in a future version of Charcoal with a function that obtains the real and/or imaginary part of a complex number:

⊞θ⁻⊟θ⊟θ≔ΣＥ²ΠＥ…θ⁺²ι§⪪ＸＩ1j∕λ⁹⁰Ｉ1jιθＩ×⁶³⁷¹⟦▷math.acosθ₂⁻²⊗θ

Don't Attempt This Online! because ATO doesn't have a new enough version of Charcoal yet. Link is to verbose version of code.

Would be only 52 bytes if AtIndex could split a complex number into real and imaginary parts:

⊞θ⁻⊟θ⊟θ≔ΣＥ²ΠＥ…θ⁻³ι§ＸＩ1j∕λ⁹⁰ιθＩ×⁶³⁷¹⟦▷math.acosθ₂⁻²⊗θ

JavaScript (ES6), 106 bytes

Same formulas (and same test case) as all other answers so far.

Expects (longitude1, longitude2, latitude1, latitude2).

with(Math)f=(p,q,r,s,R=12742)=>[R*sin(t=acos(cos(r*=x=PI/180)*cos(s*=x)*cos(x*=q-p)+sin(r)*sin(s))/2),R*t]

Try it online!

Formulas

Given the longitudes \$\lambda_1\$, \$\lambda_2\$ and the latitudes \$\phi_1\$, \$\phi_2\$:

$$t = \frac{1}{2}\arccos{(\cos{\phi_1}\cos{\phi_2}\cos{(\lambda_2-\lambda_1)}+\sin{\phi_1}\sin{\phi_2})}\\ d = R\sin t\\ s = Rt$$

where \$R=2\times 6371=12742\$

Degrees to radians conversions

In order to avoid extra parentheses, we use the *= operator for all conversions with carefully ordered operations:

cos(r *= x = PI / 180) * // define the conversion ratio x
                         // and apply the conversion to r
cos(s *= x) *            // apply the conversion to s
cos(x *= q - p) +        // x can now be trashed:
                         // this time we multiply x by q - p
sin(r) *                 // r has already been converted
sin(s)                   // s has already been converted