Weighting
The best way for image combination is a weighted mean:
If the individual weighting factors are set to wi = 1, then one has the straigh mean. If the weighting factor for a pixel is zero, then it is excluded from the coaddition (outlier rejection). The weighting scheme allows one to assign individual weights for every pixel, depending on the reliability of the informationcarried in it. Pixels with higher noise get a lower weight, hence the S/N of the stacked image is maximised. The weighting factors can directly be taken from the normalised flat. It tells which pixels are more sensitive than others, and thus carry more reliable information.
Cosmic ray rejection, satellite masking, removal of blemishes
The nice feature of this approach is that one creates one copy of the normalised skyflat for your every SCIENCE image. This WEIGHT image is then individually modified, i.e. bad or unwanted pixels are zeroed. Apart from masking bright reflections all these steps can be automatized.
Fig. 5: A WEIGHT image, with a satellite track and a few bad bixels set to zero.
Taking into account varying photometric conditions and seeing
One can even go a step further. Since your target fields rises and sets, you observe it through a different amount of airmass. Besides, the transparency of the sky can change during the observation. 0.1 mag to 0.2 mag absorption are not uncommon, in case of cirrus this can become much more. An image taken in bad photometric conditions carries more noise, which can be taken into account in the weight image. The same holds for the seeing. An image taken at a seeing of 1.5" has a S/N for faint objects that is about (2.0/1.5)2 = 1.8 better than if it was taken in 2" conditions. This too can be reflected by the individual weighting factors.
With a proper weighting scheme at hand, the S/N in a stacked image can be improved by some 10%.