Thank you JimB. My real data have usually multiple measurements for each point and I weigh the data with the variance. Mostly my models are very similar and typically only vary in initial conditions, and I expect similar error variance in the models.

Thank you Lukas, I need to study the definition a bit more but your solution looks good. I find the MultiNonlinearModelFit function very useful, and with the explicit "grad" definition you used, I hope to be able to use the "MeanPredictionBands" for a wider set of FittedModels.

Nice, this solves the issues and I do not need to give boundary constraints. But the MeshRegion of the ToBoundaryMesh has dimension 2, how do I get the volume of the region? (Strangely, the dimension of the BoundaryDiscretizeRegion, as used above, is 3. Then I can simply use the Volume function.)

Thanks, and importantly it also works if I apply BoundaryDiscretizeRegion on regInt[[I]], not only on iShell[[I]]. I want to use the function for irregular regions and must use RegionIntersection, (i.e. I cannot depend on rectangular bounds to define the intersection). Indeed, it is a workaround, but it helps me for now!

I have used Mathematica v12, on MacOS 10.14.5

OK, I added a separate paragraph with an explicit question

I want to discretize all the RegionIntersections. Eventually I want to do this for irregular regions, and estimate the volume of the intersections.