How to model and subtract nebular emission lines from stellar spectrum I have optical stellar spectra from a star forming region with spatially varying nebular background.

I want to subtract the nebular lines from the stellar spectra (which are varying in intensity from spectrum to spectrum because of the varying background).  How can I model and subtract individual nebular lines from a stellar spectrum?  Are there any known tools for this?  Any help or suggestions would be highly appreciated.
In general, you can't.  If obtaining spectra in regions where there is expected to be a spatially varying background then you either need to do long-slit spectroscopy so that you have a good measurement of the ISM contribution either side of your source, or you do integral field spectroscopy with the same idea.

The problem is that the line strengths for the emission lines you see from the ISM are integrated along the line of sight and are very dependent on the physical conditions of the ISM. This means that you cannot really assume any particular ratio of line strengths to assist in any subtraction.

In some circumstances you may be fortunate that there is a relative velocity between your stellar source and the background - in which case, separation of the nebular and stellar components may be possible. Or it may be that the stellar features are broadened by rotation (if you have the resolving power), whereas the nebular contribution is usually very narrow/unresolved.

It is difficult to provide more details without a lot more information about the observations, the target and the background.