This is the preprint of an article that is under review in the Journal of Glaciology. The abstract is as follows: Glaciological mass balance is computed from point-scale field data at a few ablation stakes that are regressed as a function of elevation, and averaged over the area-elevation distribution of the glacier. This method is contingent on a tight control of elevation on local ablation. On debris-covered glaciers, systematic and random spatial variations of debris thickness modify ablation rates. A method that takes into account the debris-thickness variability in extrapolating point-scale ablation data may be more accurate on these glaciers. We propose and test a methGlaciological mass balance is computed from point-scale field data at a few ablation stakes that are regressed as a function of elevation, and averaged over the area-elevation distribution of the glacier. This method is contingent on a tight control of elevation on local ablation. On debris-covered glaciers, systematic and random spatial variations of debris thickness modify ablation rates. A method that takes into account the debris-thickness variability in extrapolating point-scale ablation data may be more accurate on these glaciers. We propose and test a method where stake data are interpolated as a function of debris-thickness alone, and averaged over the observed debris-thickness distribution at different parts of the glacier. We apply this method to compute sub-debris ablation rate on Satopanth Glacier (Central Himalaya) utilising about a thousand ablation measurements at a network of up to 56 stakes during 2015--2017. We compare our results with that from the standard glaciological method. The uncertainties in both the estimates due to the corresponding uncertainties in measurement of ablation and debris-thickness distribution, and that due to interpolation procedures are estimated using Monte Carlo methods. Possible biases due to finite number of stakes used are investigated, and net specific balance of Satopanth glacier is computed.od where stake data are interpolated as a function of debris-thickness alone, and averaged over the observed debris-thickness distribution at different parts of the glacier. We apply this method to compute sub-debris ablation rate on Satopanth Glacier (Central Himalaya) utilising about a thousand ablation measurements at a network of up to 56 stakes during 2015--2017. We compare our results with that from the standard glaciological method. The uncertainties in both the estimates due to the corresponding uncertainties in measurement of ablation and debris-thickness distribution, and that due to interpolation procedures are estimated using Monte Carlo methods. Possible biases due to finite number of stakes used are investigated, and net specific balance of Satopanth glacier is computed.