OK, here is what I would have done.
Key things I did, were use the color propagation method of highlight recovery, like others did, and then use filmic tone curves and the CL chromacity according to lightness curve to fade the highlights to white in an organic way. Other things I did to benefit this high dynamic range scene, were, Dynamic Range compression, Graduated filter lightening the shadows, Shadows highlights with heavy shadows lifting, a boost of local contrast using wavelets, especially some of the lower frequencies to keep things from getting too salt and pepper-ish, As I saw lots of interesting color detail yet there was excessive chroma noise, I applied conservative rgb noise reduction and added some of the color back with low frequency chroma wavelet boosting.