If we assume that light travels in straight lines we could easily recognize the real shape of the earth and calculate its size by using a theodolite to measure the drop of a target or the horizon. On a flat earth the drop is zero on the globe it is about d^{2} / 2R, where d = distance to target and R = radius of the earth.
But we have an atmosphere which refracts light up or down. What a theodolite sees is not the real geometry of the earth, but an up or down refracted image. Whatever the shape of the earth is, we will always see a distorted image of the earth due to refraction. So how can we measure how much curvature is due to the shape of the earth and how much is due to refraction?
If we can measure refraction somehow, we can correct the apparent curvature and get the real geometry of the earth.
Note: there are devices that can measure refraction directly, e.g. by using 2 lasers of different wavelengths. So we know the earth is a globe not at last from such measurements. [1]
Refraction depends on the atmospheric condition between observer and target. We can approximate the light rays between observer and target through an arc with a radius of curvature r.
We can use the method of simultaneous reciprocal zenith angle measurements using 2 theodolites to measure the apparent drop the theodolites see each other. The drop is measured as the angle from the vertical to the target. The angle is called the Zenith Angle. The sum of the measured zenith angles depends on the shape of the earth and the radius of curvature of the light ray. So if we know the shape and size of the earth, we can calculate the refraction from the zenith angle measurements.
Now we pretend not to know the shape of the earth. But we can calculate the curvature of the light ray for both models anyway and then use a physical model of the atmosphere to calculate the atmospheric parameters that result in this curvature. Then we can measure the real parameters like the temperature, pressure and temperature gradient and show for which earth model the real values can match the calculations.
If the real conditions can't produce the curvature of the light that is calculated for one or the other earth model, than that model is false.
So lets see how this works.
Refraction calculations commonly assume the earth is a globe. This is expressed in the refraction coefficient k which is defined as:
(1) 

per definition  
where^{'} 

Because refraction itself does not depend on the shape of the earth, we can derive all calculations without assuming the shape of the earth.
We first need the connetion between the curvature of a light ray and the corresponding refraction angle. Then we can derive equations for the refraction angle from the zenith angles for globe and flat earth model. Lets start:
The light ray curvature is not dependent on the earth model and can be calculated from the refraction angle as follows.
The magenta line is the direct line of sight from the observer to a target. The orange arc is the corresponding refracted light ray with radius of curvature r. The refraction angle is ρ (greek r, spelled rho) and the distance to the target is d.
Although the sketch shows a globe model, the light ray geometry has nothing to do with the shape of the earth. The radius of the earth R does not appear in the equations below.
We have a right angle triangle r, d/2, b with an angle ρ. We can use trigonometry to calculate the sine of ρ:
(2) 
The curvature κ of the bent light rays is very small so the radius of curvature r = 1 / κ is very big. We have very small refraction angles ρ if d is much smaller than r, which is the case in survey. For small angles ρ we can approximate sin(ρ) ≈ ρ:
(3) 
The curvature of the light ray is the inverse of the radius of curvature. So we can solve (3) for 1/r:
(4) 
 
where^{'} 

Now we derive the equation for the refraction angle
Assuming we have measured both zenith angles and assuming the radius of the earth is R we can get the refraction angle
(5) 
Rearranging for
(6) 
The angle
(7)  
where^{'} 

Inserting this into (6) and solving for the refraction angle
(8) 
 
where^{'} 

We can now insert the refraction angle into the equation (4) for the curvature of a light ray:
(9) 
This can be simplified to:
(10) 
 
where^{'} 

This is the equation to calculate the curvature of a light ray from the simultaneously measured zenith angles for the globe model.
We can get the commonly used refraction coefficient k as follows:
(11) 

Note: if the light ray is straight, the right term above is equal −1 so the refraction coefficient for unrefracted light is k = 0.
We can use equation (11) for the flat earth model with a little correction. Because the verticals at the 2 locations on the flat earth are parallel, the tilt angle
(12) 
 
where^{'} 

Similar as above we can now calculate the equation for the curvature of a light ray from the simultaneously measured zenith angles for the flat earth model:
(13) 
 
where^{'} 

This is alomost the same equation as for the globe model (10), but the 1 / R term does not appear.
We can even calculate a refraction coefficient for the flat earth with the same definition
(14) 

which is the same as the refraction coefficient for the globe model (11) minus 1. So whenever we have a refraction coefficient for the globe, e.g. in the Advanced Earth Curvature Calculator, we can calculate the corresponding flat earth refraction coefficient as:
(15) 

Note: in reality the sum of the zenith angles is almost always measured greater than 180°. This means the curvature of the light ray on the flat earth model is negative, so the atmospheric conditions have to be such, that light gets bent upwards. This means that the air density has to increase with increasing altitude, which is not the case in reality.
But lets investigate the atmospheric conditions neccessary for a given zenith angle observation for the globe and the flat earth:
The curvature of the light in the atmosphere is due to a gradient in the refractive index. The refractive index depends on air density, which is dependent via the ideal gas law on pressure and temperature, on the wavelenght of light and to a very small extent on humidity and CO2 concentration.
Because refraction does fluctuate in practice all the time, for practical calculations humidity and CO2 concentration can be ignored and the wavelength of green light is used. Try my Calculator for Refractivity based on Ciddor Equation and play with the sliders to see how much the influence of each parameter on the refractive index is.
Knowing the connection between atmospheric conditions and the refractive index (Ciddor Equation), using calculus we can derive the curvature of a light ray passing through an atmosphere with a certain density gradient, see Deriving Equations for Atmospheric Refraction.
The following equation gives the local curvature of a light ray depending on the atmospheric conditions at the location. If the conditions are about the same from the observer to the target, we can use this curvature for the whole path. In geodesy the inclination angle of a light ray is always near zero, so the influence of this angle (factor cos(≈0°) = 1) can be ignored.
(16) 
 
where^{'} 

Note that this equation does not depend on the shape of the earth. Refraction is independent of the shape of the earth. So we can use this equation for the globe and flat earth model.
I have shown how we can calculate the curvature of light from measurements of the zenith angles for the globe and flat earth model.
If we solve the equation (16) above for the temperature gradient, we can then insert the calculated light ray curvatures to calculate the temperature gradient neccessary to bend the light ray accordingly:
(17) 
 
where^{'} 

This equation is independent of the shape of the earth as well.
To calculate the globe an flat earth results the Calculator below was used. The elevation (elev) is the mean elevation of the 2 stations. It's not the observer height above the ground. The Barometric Calculator was used to calculate pressure and temperature for this mean elevation. The Baromentric Calculator assumes the International Standard Atmosphere.
Note: The observer height is very important to get out of the strong refraction layer, see Refraction Coefficient as a Function of Altitude. For observations over great distances not over a flat landscape or water the line of sight is commonly for the most part way above the ground layer, so we can expect standard refraction, as the values in the table below confirm. On short distance measurements with observer heights around 2 m, the light rays travel in most cases the whole distance in the ground layer, where refraction can be considerably different than standard in both direction, which is confirmed by the Problematic Observations.
Measurements  Globe Model  Flat Earth Model  

Src    dist. [m]  elev. [m]  ρ  r [km]  k  dT/dh [°C/km]  ρ  r [km]  k  dT/dh [°C/km] 
NGS  91°08'09"  89°07'44"  32,139  2273  43.8"  75,745  0.084  −18.0  7'56.5"  −6956  −0.916  −211 
NGS  91°03'44"  89°15'34"  39,476  2188  1'00.0"  67,821  0.094  −16.3  9'39.0"  −7032  −0.906  −208 
NGS  90°40'36"  89°56'32"  76,957  2308  2'11.8"  60,235  0.106  −13.8  18'34.0"  −7125  −0.894  −208 
NGS  90°32'21"  90°12'12"  92,882  2421  2'47.1"  57,342  0.111  −12.6  22'16.5"  −7167  −0.889  −208 
NGS  90°49'49"  89°45'41"  74,362  2209  2'18.8"  55,270  0.115  −12.1  17'45.0"  −7201  −0.885  −205 
NGS  90°50'01"  89°35'13"  52,518  2269  1'33.2"  58,146  0.110  −13.1  12'37.0"  −7155  −0.890  −206 
NGS  91°09'00"  89°22'35"  66,128  2043  2'03.0"  55,462  0.115  −12.5  15'47.5"  −7198  −0.885  −202 
JK1  90°00'33"  90°00'34"  2228  0.84  2.6"  89,324  0.071  −22.6  33.5"  −6860  −0.929  −186 
JK1  90°00'29"  90°00'27"  2228  0.84  8.1"  28,468  0.224  2.16  28.0"  −8208  −0.776  −161 
LS  88°57'58.2"  91°09'55"  17,824  387  51.9"  35,397  0.180  −4.13  3'56.6"  −7769  −0.820  −172 
This observations were done over great distances with a line of sight over 30 m above the ground, so the temperature gradient was not influenced by the ground.
For the flat earth model we get temperature gradients of −161°C/km to −211°C/km. Such steep gradients are only possible for light rays traveling very low along a hot surface with cool air above. The apparent globular shape of the earth demands such a gradient in any height over the ground on the flat earth. This is physically impossible because the temperature would reach absolute zero below 2 km altitude.
On the globe model however the calculated temperature gradients are as expected. According to a realistic atmospheric model the gradients will fade to about −6.5°C/km of the standard atmosphere in the lower 30 m above the ground and maintain this gradient up to about 11 km, see Refraction Coefficient as a Function of Altitude. Due to decreasing air density with increasing altitude, refraction decreases with altitude even with a constant temperature gradient. The range of refraction of k = 0.07..0.2 with a mean value of 0.14 corresponds to standard refraction.
If we measure above 30 m ground level, we get consistent refraction with not much variation, see Refraction Coefficient as a Function of Altitude.
This measurements show that the earth cannot be flat. The temperature gradients neccessary to produce this observations on the flat earth are physically impossible. The calculated temperature gradients for the globe on the other hand are in the neighborhood of the gradient −6.5°C/km of the standard atmosphere.
If the observer measures over relatively short distances from a height of about 2 m, the light travels in the ground layer of the atmosphere, where the ground exchanges heat with the air above, creating steep temperature gradients and hence strong refraction.
If the ground is cooler than the air, e.g. in the evening when the ground cools down faster than the air above, the gradient is always strong enough to bend light along the surface for hundreds of km. So it is no surprise, that on the globe earth we can see lasers from behind the curvature in any distance. The nearer to the ground the observer and laser are, the stronger refraction is.
This ground effect is very well supported by the measurements in the following table, where all observations were done over short distances (dist) from below 30 m observer height (ht). This measurements are inconclusive to tell the shape of the earth, although the flat earth temperature gradients are much more pronounced.
Measurements  Globe Model  Flat Earth Model  

Src    ht. [m]  dist. [m]  elev. [m]  ρ  k  dT/dh [°C/km]  ρ  k  dT/dh [°C/km] 
JK2b  90°00'11"  90°01'00"  1.7  1538  16.5  10.60"  −0.426  −104  35.50"  −1.426  −267 
JK2c  89°40'31"  90°19'58"  1.6  1100  114  3.31"  0.186  −3.8  14.50"  −0.814  −168 
JK2c  90°22'18"  89°37'51"  1.6  150  110  2.07"  −0.853  −174  4.50"  −1.853  −339 
JK2c  89°43'56"  90°16'18"  1.6  275  110  2.55"  −0.572  −128  7.00"  −1.572  −293 
JK2d  90°14'24"  89°45'43"  1.6  143  11.8  1.19"  −0.512  −118  3.50"  −1.512  −281 
JK2d  89°55'25"  90°04'41"  1.6  129  12  0.91"  −0.437  −105  3.00"  −1.437  −269 
JK2e  90°04'32"  89°55'33"  1.6  61  5.6  1.51"  −1.536  −285  2.50"  −2.536  −448 
JK2e  87°57'53"  92°02'12"  1.6  61  6.7  1.52"  −1.538  −285  2.50"  −2.538  −448 
JK2f  90°21'16"  89°38'48"  1.6  150  28.4  0.43"  0.176  −5.33  2.00"  −0.824  −169 
JK2f  90°35'41"  89°25'27"  1.8  1000  23.7  17.81"  −1.10  −234  34.00"  −2.10  −377 
JK2f  90°59'21"  89°00'47"  1.6  250  25.7  0.05"  0.012  −32.3  4.00"  −0.988  −196 
JK2f  90°33'15"  89°27'04"  1.6  500  21.1  1.40"  −0.173  −62.6  9.50"  −1.17  −226 
JK2f  90°27'37"  89°32'47"  1.7  600  21.1  2.29"  −0.235  −72.7  12.00"  −1.24  −236 
JK2f  89°59'04"  90°01'03"  1.6  100  18.8  1.88"  −1.16  −224  3.50"  −2.16  −387 
JK2i  90°00'22"  90°00'04"  1.6  1538  16.5  11.9"  0.478  43.7  13.0"  −0.522  −119 
JK2i  90°00'09"  89°59'54"  1.6  1538  16.5  23.4"  0.940  119  1.5"  −0.060  −441 
Note: Because the refraction error increases with the square of the distance, for short enough distances the error l is small enough to be not a problem.
(18) 
 
where^{'} 

We can calculate the maximal distance to assure a certain accuracy even with high refraction. Measuring multiple times at different times a day and average over the measurements can give very good results. To do so we solve (18) for d and set l to the maximal error we accept.
The calculator is preset for the first observation JK1 above. You can inspect the Calculator Code below.
ID  Date  Project  Description 

a  1977  NGS  McDonald Observatory  Long distance (32..93 km) observations at about 2000 m elevation, line of sight above 30 m 
NGS Zenith Angle Measurements from the McDonald Observatory (PDF)
Survey of the McDonald Observatory Radial Line Scheme by Relative Lateration Techniques
NOAA Technical Report (PDF), June 1978
The NOS/NGS performed a special survey in the vicinity of the University of Texas McDonald Observatory. This was the initial phase of an extensive geodeticgeophysical study to detect any motions of the observatory relative to prominent topographic features within a region extending as far as 100 km from the observatory.
Very Long Line EDM & Reciprocal Zenith Angle Observations – A Measure Of The Plumbline Tilt
Because that report included their reciprocal zenith angles (they call them zenith distances), I chose to include that set (see JK2a) and I personally knew Charlie Glover who performed that work. Here is a blog entry I made to share that information. ~ Jesse Kozlowski
Flat Earth Proof? – Perfectly Flat Level Lake by Jesse Kozlowski
Detailed documentation of the measurements and all data of the observation
ID  Date  Project  Description 

a  1977  NGS  McDonald Observatory  Long distance (32..93 km) observations at about 2000 m elevation, line of sight above 30 m 
b  NOV 09 2015  RTE 206 MP 27 & 28  d = 1500 m, observer height 1.6 m 
c  MARCH 08 2016  Clarksville EDM CBL  d = 275..1000 m, observer height 1.6 m 
d  MARCH 30 2016  HAINESPORT "L" SHAPE  d = 130..142 m, observer height 1.6 m 
e  APRIL 20 2016  200 FOOT "L" SHAPE  d = 60 m, observer height 1.6 m 
f  OCT 07 2016  Folsom EDM CBL S8  d = 100..1000 m, observer height 1.6 m 
g  OCT 20 2016  Union Lake GNSS T2  d = 760 m, observer height 1.4 m 
h  NOV 01 2016  Cooper River Lake  d = 2200 m, observer height 1.6 m 
i  NOV 617 2016  RTE 206 MP 27&28 LEVEL LINE  d = 1500 m, observer height 1.6 m 
var Model = { R: 6371000, z1: 90.00916667, z2: 90.00944444, d: 2228.4, P: 1013.25, // mbar T: 15, // °C rho_gl: 0, rho_fe: 0, r_gl: 0, r_fe: 0, kappa_gl: 0, kappa_fe: 0, k_gl: 0, k_fe: 0, dTdh_gl: 0, dTdh_fe: 0, Update: function() { // globe model var TK = this.T + 273.15; var a = (180  (this.z1 + this.z2)) * Math.PI / 180; this.rho_gl = this.d / 2 / this.R + 0.5 * a; this.kappa_gl = 2 * this.rho_gl / this.d; this.r_gl = 1 / this.kappa_gl; this.k_gl = this.kappa_gl * this.R; this.dTdh_gl = 12666 * this.kappa_gl * TK*TK / this.P  0.0343; // flat earth model this.rho_fe = 0.5 * a; this.kappa_fe = a / this.d; this.r_fe = 1 / this.kappa_fe; this.k_fe = this.kappa_fe * this.R; this.dTdh_fe = 12666 * this.kappa_fe * TK*TK / this.P  0.0343; ControlPanels.Update(); }, }; ControlPanels.NewPanel( { Name: 'InputPanel', ModelRef: 'Model', NCols: 2, OnModelChange: function(field) { Model.Update(field) }, Format: 'std', Digits: 5, ReadOnly: false, PanelFormat: 'InputNormalWidth' } ).AddHeader( { Text: 'Enter Zenith Angles, Distance, Pressure and Temperature', ColSpan: 4, } ).AddTextField( { Name: 'z1', Label: 'ζ<sub>1</sub>', Format: 'dms', ConvToModelFunc: function(s){ return NumFormatter.DmsStrToNum(s); }, Digits: 5, } ).AddTextField( { Name: 'z2', Label: 'ζ<sub>2</sub>', Format: 'dms', ConvToModelFunc: function(s){ return NumFormatter.DmsStrToNum(s); }, Digits: 5, } ).AddTextField( { Name: 'd', Label: 'Distance', Units: 'm', } ).AddTextField( { Name: 'R', Label: 'Radius Earth', Units: 'km', Mult: 1000, } ).AddTextField( { Name: 'P', Label: 'Pressure', Units: 'mbar', Digits: 6, } ).AddTextField( { Name: 'T', Label: 'Temperature', Units: '°C', } ).Render(); ControlPanels.NewPanel( { Name: 'OutputPanel', ModelRef: 'Model', NCols: 2, OnModelChange: function(field) { Model.Update(field) }, Format: 'std', Digits: 5, ReadOnly: true, PanelFormat: 'InputNormalWidth' } ).AddHeader( { Text: 'Results', } ).AddHeader( { Text: 'Globe', } ).AddHeader( { Text: '', } ).AddHeader( { Text: 'Flat Earth', } ).AddTextField( { Name: 'rho_gl', Label: 'Refr.Angle ρ', Mult: Math.PI / 180, Format: 'dms', Digits: 6, } ).AddTextField( { Name: 'rho_fe', Label: 'ρ', Mult: Math.PI / 180, Format: 'dms', Digits: 6, } ).AddTextField( { Name: 'kappa_gl', Label: 'Light Curve κ=1/r', Format: 'sci', Units: '/m', } ).AddTextField( { Name: 'kappa_fe', Label: 'κ', Format: 'sci', Units: '/m', } ).AddTextField( { Name: 'r_gl', Label: 'Light Curve r', Units: 'km', Mult: 1000, } ).AddTextField( { Name: 'r_fe', Label: 'r', Units: 'km', Mult: 1000, } ).AddTextField( { Name: 'k_gl', Label: 'Refr.Coeff k', } ).AddTextField( { Name: 'k_fe', Label: 'k', } ).AddTextField( { Name: 'dTdh_gl', Label: 'Temp.Gradient dT/dh', Units: '°C/km', Mult: 0.001, } ).AddTextField( { Name: 'dTdh_fe', Label: 'dT/dh', Units: '°C/km', Mult: 0.001, } ).Render(); xOnLoad( function() { Model.Update(); } );