Background: In planetary astronomy, the only method we have to estimate the age of a surface of a solid object in the solar system (other than Earth) is to identify craters. We then compare these to an established chronology that links the number density of craters with diameters $D \ge 1$ km (often written as $N(1)$) with absolute ages as established from Apollo and Luna sample return missions.

The function linking age $T$ with $N(1)$ has historically been fit to the form:

$N(1) = \alpha (\exp(-\beta T)-1)+\gamma T$

where $\alpha$, $\beta$, and $\gamma$ are parameters of the fit. Qualitatively, this expresses the idea that early on, there was an exponential decrease in cratering, and after a certain time, there was a linear rate up to the present day.

This fit is constrained by around 15 points at most, it's completely unconstrained (other than the maximum mass of an object) for ages older than about 3.95 billion years, and there are no data points between about 1 and 3 billion years ago.

Question: I'm re-doing some of this older work as part of my research now. I'm using established radiometric ages and re-doing the crater counts to fit a new function. My issue is that, from my new data, this version of the fit function "looks" bad. A research group in 2007 suggested more of a quadratic decline in the recent past so revised the function as:

$N(1) = \alpha (\exp(-\beta T)-1)+\gamma T^2 + \delta T$

Adding the quadratic term greatly increases how good the fit "looks," but with only 11 points, I'm worried that I could be getting into an area where I'm just getting an improved fit because I'm adding more free parameters. And, probably, a reviewer is going to want something more quantitative, especially because the paper that everyone uses for this has the first version.

Note that the data range is roughly $0 < T < 4.5$ and $10^{-6} < N(1) < 10^0$, and while they're somewhat well-behaved, it's non-linear (as the function would suggest).

I've calculated the reduced $\chi^2$ of each version, and it's 2.8 versus 1.7, respectively. But, this is a highly non-linear function, and based on some reading I've done, the $\chi^2$ may not be a meaningful metric for determining how good the fit really is. Someone suggested that I do an incomplete gamma function test to determine if the $\chi^2$ is meaningful, but the results of that are 0.00018 versus 0.017 ... and I have no idea what those mean other than I had a thought that the larger (if still minuscule) number was better indicating the quadratic was more meaningful.

So ... what's a (or several) good way to determine whether one fit function is statistically better than another with this kind of data? $\chi^2$? Or something else?