Thermal membrane distillation (MD) is an emerging technology to desalinate high-salinity wastewaters, including shale gas produced water to reduce the corresponding water footprint of fracturing operations. In this work, we introduce a rigorous optimization model with energy recovery for the synthesis of multistage direct contact membrane distillation (DCMD) system. The mathematical model (implemented in GAMS software) is formulated via generalized disjunctive programming (GDP) and mixed-integer nonlinear programming (MINLP). To maximize the total amount of water recovered, the outflow brine is fixed close to salt saturation conditions (300 g.kg(-1) water) approaching zero liquid discharge (ZLD). A sensitivity analysis is performed to evaluate the system's behavior under different uncertainty sources such as the heat source availability and inlet salinity conditions. The results emphasize the applicability of this promising technology, especially with low steam cost or waste heat, and reveal variable costs and system configurations depending on inlet conditions. For a produced water salinity ranging from 150 g.kg(-1) water to 250 g.kg(-1) water based on Marcellus play, an optimal treating cost are between 11.5 and 4.4 US$ m(-3) is obtained when using low-cost steam. This cost can decrease to 2.8 US$ m(-3) when waste heat from shale gas operations is used.