Orthogonal coordinate systems represent an attractive alternative for the formulation of two-dimensional composition space equations for partially premixed combustion: They avoid the need of closure models for a cross scalar dissipation rate and allow for a direct recovery of the corresponding flamelet equations in the asymptotic limits of non-premixed and premixed combustion. Despite these remarkable features, this kind of coordinate system still presents some important unsolved issues, which have limited their application so far. These difficulties are mainly associated with i) the lack of an appropriate formalism for the definition of two variables having orthogonal gradients in the entire flame domain and ii) the absence of corresponding closure models for the gradients of these variables in two-dimensional composition space. In the present work, it is shown how a Lagrangian interpretation of the flamelet derivative allows solving both problems. More specifically, after the mixture fraction, Z, is adopted as first coordinate, the proposed approach allows to derive i) two-dimensional composition space equations for all reactive scalars, ii) a transport equation for a modified reaction progress variable, φ, satisfying the desired orthogonality condition, ∇Z·∇φ=0, and iii) two-dimensional composition space equations for gZ=|∇Z| and gφ=|∇φ|. The obtained set of two-dimensional equations in orthogonal composition space is general and it can describe different flames of interest. In order to illustrate the capabilities of the formulation, the equations are then specialized for planar flames with unity Lewis number, obtaining in this way a solvable set of composition space equations in orthogonal coordinates. After an appropriate numerical approach is introduced, the resulting 2D equations are solved to analyze the interaction between premixed flamelets with a strain rate prescribed along the Z-dimension controlling the interaction. Both flame structures and budgets of the scalar gradient equations are studied and the results provide new insights into the physics of scalar gradients in two-dimensional composition space. Finally, conceivable coupling strategies of the present formulation with CFD codes for the simulation of turbulent flames are discussed.