Here we describe the coupling of the chemical speciation software ORCHESTRA with the parameter estimation software PEST. This combination enables the computation of optimised model parameters from experimental data for the ion binding models implemented in ORCHESTRA. For testing this flexible tool, the NICA-Donnan model parameters for proton-, Cd- and Zn-binding to Laurentian fulvic acid were optimised. The extensive description of the method implementation and the examples provided facilitate the use of this tool by students and researchers. Three procedures were compared which derive the proton binding parameters, differing in the way they constrain the model parameters and in the implementation of the electrostatic Donnan model. Although the different procedures resulted in significantly different sets of model parameters, the experimental data fit obtained was of similar quality. The choice of the relation between the Donnan volume and the ionic strength appears to have a strong influence on the derived set of optimal model parameters, especially on the values of the protonation constants, as well as on the Donnan potential and Donnan volume. Optimised results are discussed in terms of their physico-chemical plausibility. Coherent sets of NICA- Donnan parameters were derived for Cd and Zn binding to Laurentian fulvic acid.