2D Profile, multiple aquifers

We have developed a 2D profile simulation of layered aquifers within inflows/outflows controlled by prescribed heads.  Although our interest is in the phreatic surface, the model forces confined conditions regardless of the configuration.  We have converted to unsaturated flow, still the results show full saturation of the domain, suggesting confined conditions. Short of converting to 3D, is there a 2D solution?

Unless doing an unconfined simulation, a 2D vertical cross-section / profile model is always calculated as saturated. So far your solution is fine. Knowing why the unsaturated simulation leads to saturated conditions is difficult. One possibility is the coordinate reference. For cross-sectional models, the reference of the vertical y coordinate is important. Internally, FEFLOW always uses the local coordinate system, so you should make sure that there is no offset between the local and global systems.