Adding band to existing GeoTiff using GDAL?

The solution, if the driver supports it, is to call GDALOpen() with GA_Update access then use GDALAddBand or GDALDataset::AddBand. However, the geotiff driver doesn't support AddBand.

To expand on Luke's answer and provide a concrete example in Python, here's a snippet that adds an alpha band to a source raster and saves it as a PNG.

from osgeo import gdal

src_ds = gdal.OpenShared(input_path)
mask_ds = gdal.OpenShared(mask_path)
mask = mask_ds.GetRasterBand(1).ReadAsArray()

tmp_ds = gdal.GetDriverByName('MEM').CreateCopy('', src_ds, 0)
tmp_ds.AddBand()
tmp_ds.GetRasterBand(4).WriteArray(mask)

dst_ds = gdal.GetDriverByName('PNG').CreateCopy(output_path, tmp_ds, 0)
del dst_ds

I used the MEM driver instead of VRT since the latter does not support WriteRaster() and WriteArray() (error "Writing through VRTSourcedRasterBand is not supported."). Using the vrt driver might still be possible through some other methods, I suppose.