We use a modified form of the Burke et. al (2020) bottom hole temperature (BHT) method on a county-scale basis. The concept is that every BHT is too low because of unknown circulation and shutin times, as well as unrecorded mud temperature after circulation. Our modification is to anchor the gradients at the 10-yr average surface temperature, set up a sweep of gradient lines on the high side of the BHTs, then visually determine the highest geothermal gradient (G) supported by the data. This G forms a lower bound on the geothermal gradient to the deepest support points, the gradient could be better (hotter) but not worse (cooler). Below is an example for Columbia County, Arkansas, where the data supports a gradient of > 40 C/km at 2500 m. 

We suspect that the gradient is about the same for the 3500 m zone, but have insufficient number of BHTs to prove it. This is a statistical approach and it takes a large number of BHTs to have a few that, through serendipity, approach the true geothermal gradient.

You might notice the very high gradients implied shallower in the plot. We see in all areas that G decreases with depth: in other words, 'the geothermal gradient' does not exist, there are many Gs at any given location. Geothermal gradient is a multivalued (nonlinear) function of depth. Even so, a linear geothermal gradient is a useful approximation for estimating reservoir temperature for a specific reservoir depth.