This manuscript outlines a new optimized method for constraining Multiple Diffusion Domain (MDD) model parameters from step-wise degassing experiments. The authors argue that the most commonly used method (Lovera et al. 1997) could be yielding too low temperature predictions due to the assumption that each domain loses no more than 60% of its gas during the lower temperature heating steps. The user can introduce bias when fitting a regression to the lower temperature steps when quantifying the activation energy, which in turn can yield lower temperature predicted thermal-histories. The open-source software presented here (MDD Tool Kit) does not need user inputs for regression fits, and so does not rely on any assumptions about domain size distribution. To test the software, the authors then modelled a previously published dataset (Wong et al., 2023) using the MDD Tool Kit. The resulting thermal histories are consistent with the geological and independent constraints, but are consistently higher temperature to those predicted by Wong et al. 2023.

I find that this manuscript is a good fit for GChron. It is well written with clear figures and tables. The authors lay out problems they see with the current most used method, and describe the model well. I have a few general and specific relatively minor comments. 

General points:

1. The authors mention the criticism of the MDD model – it would be valuable to include a brief section discussing both the criticisms and supporting evidence of the model.
2. MDD theory has many potential problems – e.g., the assumption that phases remain constant during heating experiments. While it is out of the scope of this work to show evidence for or against various assumptions, the manuscript would benefit from a section including more detail on what the different assumptions actually are and to include references to previous studies on how these have been tested.
3. Which choice of misfit method do you think is best to use? How would you advise a user of this tool kit to evaluate which is the best method to use for their data? In some cases, the results can be pretty significantly different (e.g., the high temperature portion of GR27’s cooling history).
4. There are a lot of variables in this study, a table defining them could be helpful for readers to refer back to.

Specific points:

1. Figure 1 A and B the x axis ticks are not the same in the two plots.
2. Figures 3 and 5 are not directly referenced in the text.
3. The x axis for Figure 4 cooling histories are not the same – it might be quicker for readers to evaluate the cooling history vs stratigraphic depth if these were the same.
4. The synthetic data results and discussion section is a little unclear to me. You say the MDD tool kit did not perform as well for Experiment B, but it looks (in table 3) that experiment B’s results are closer to the setup values.
5. I think the final column for the supplement tables (the label for whether a step should be included in the optimization) is not correct – legend says 2 = not used, but the table the value is 0.
6. I like that you include a parameter range larger than realistic to avoid any potential bias.
7. Did you choose 30,000 iterations when modelling the data from Wong et al. (2023) to allow for complete convergence – Figure 2 suggests this happens earlier at ~18,000 iterations. How would you advise users to choose the number of iterations?
8. The cooling history for GR1 is initially pretty high temperature, especially when you compare it to the pre-existing BtAr date. The cooling history from Wong et al., 2023 would suggest that the sample experienced relatively monotonic cooling since BtAr closure. Is this realistic?

The authors provide a useful update to the semi-manual multi-domain modelling (MDD) approach originally set out in Lovera et al. (1997) and updated in subsequent papers by the same Lovera and colleagues. The manuscript focusses on optimization of the MDD model parameters and does appear to improve reproducibility of the output thermal histories while also revealing higher temperature paths by 60-75°C. The elimination of user-defined E_a and user choice of points for regression is to be applauded, it was always a significant source of potential user bias.

However, there are some issues that require further discussion and an issue with the choice of natural world example for comparison.

Firstly, there continues to be debate about the physical validity of the multi-domain modelling approach and that needs more discussion in this article. Recent challenges in the literature appear to be missing from the references such as Spikings and Popov 2021 (and references therein) which deserve mention. I will also briefly discuss the issue of 39Ar recoil below.

Spikings, R.A. and Popov, D.V. 2021. Thermochronology of alkali feldspar and muscovite at t > 150 °c using the 40Ar/39Ar method: A review, Minerals, 11, 2021 Article number 1025

Secondly, the study chosen (Wong et al. 2023, figure 4) does not demonstrate better alignment of the new python MDD code with the data Wong et al. (2023) based on biotite Ar-Ar, ZFT, AFT or Apatite (U/Th)/He. Indeed the Wong et al. thermal histories for sample GR-8 in figure 4 show a better fit for the Lovera et al. (1977) MDD model than the new approach. There is also little in the GR-1 data that indicates the new approach gives a better fit to the other constraints on natural samples. It is likely that the Ar-Ar, FT and (U,Th)/He data have flaws but the comparison with a natural example does not confirm a more realistic fit for the new model.

Line by line comments:

Line 71: Parsons et al. (1999) also discuss the physical form of alkali feldspar diffusion domains in detail.

Line 84: is any account taken of the potential for release of Ar from fluid inclusions at low temperature?

Line 95: Note that some iterations of the MDD model used variable E_a see (Lovera, O.M., Grove, M. and Harrison, T.M. 2002. Systematic analysis of K-feldspar Ar-40/Ar-39 step heating results II: Relevance of laboratory argon diffusion properties to nature. Geochimica Et Cosmochimica Acta, 66, 1237-1255.)

Line 141: add X_best to this line in explanation of the equation on line 143.

Line 145: Sentence contains typos and check square of curved bracket in the same line.

Line 146: Please explain why 0.7 is used.

Line 162: Can the validity of using 0.25 difference as the minimum difference to determine the requirement for an additional domain be tested? For example, what happens when 0.1 or 0.5 are used?

Line 171: Assuming the largest domain is similar to the size of the physical mineral grain the new approach (similar to earlier MDD models) indicates the smallest domains are commonly less than 1 micron. Given the known mineral and textural features of alkali feldspar of similar dimensions (e.g. Parson et al. 1999), please comment on the potential for 39Ar recoil affecting the estimates of the lower end of the output thermal history using the new approach.

Line 185: Authors are commended for including this physical info on feldspar melting, note that some of the early MDD modelling did not exclude such high temperature steps and thus direct comparison with the new model will be challenging.

Line 254: The authors present the MDD tool kit as an unbiased approach but at this point they leave the reader or MDD model practitioner with two ways to optimise the model, stating that the ‘true’ E_a lies between that defined by χ­² or %_frac . Please add discussion to aid the reader or practitioner on how to discriminate and choose the best estimate. Without this discussion the way is open for the practitioner to choose the value that gives the thermal history they prefer, rather than the one best fit by the model.

Line 284: As above in the initial comments, I agree the new MDD tool kit is producing more reproducible and consistent outputs. However, the new thermal histories do not fit the few independent temperature/time constraints shown in figure 4 better than the earlier modelling. Can the authors discuss and point to where the fit is better or the new thermal history is significantly better fit to the natural constraints of the Grayback fault area.

On balance the MD tool kit is an improvement on earlier approach, giving more consistent and reproducible outputs. While the comparison using a set of natural examples is not a clear endorsement of the new approach the authors should acknowledge that more work on natural examples are required to fully test the approach.