The authors consider two incompressible fluids with constant densities \(\rho _{-}\), \(\rho _{+}\) and equal viscosity \(\mu \), separated by a connected curve \(z^{o}=(z_{1}^{o},z_{2}^{o})\) inside a 2D porous medium with constant permeability \(\kappa \), and submitted to gravity \(-g(0,1)\). They write the incompressible porous media system as \(\partial _{t}\rho =\nabla \cdot (\rho v)\), \(\nabla \cdot v=0\), \(\frac{\mu }{\kappa }v=-\nabla p-\rho g(0,1)\), where \(v(t,x)\) is the velocity field and \(p(t,x)\) the pressure. The paper focuses on two unstable solutions to this Muskat problem where the two fluids occupy initial domains \(\Omega _{\pm }^{o}\), which are separated either by a closed chord-arc curve, in the case of so-called bubble interfaces, or by an open chord-arc curve which cannot be parametrized as a graph, in the case of turned interfaces. In the bubble interface case, the authors prove that for every closed chord-arc curve \(z^{o}\in H^{6}(\mathbb{T };\mathbb{R}^{2})\) there exist infinitely many mixing solutions to the Muskat problem starting from \(\rho ^{o}(x)=\rho _{\pm }=\pm 1\), \(x\in \Omega _{\pm }^{o}\). In the turned interface case, the authors prove that for every open chord-arc curve, either \(x_{1}\)-periodic \(z^{o}-(\alpha ,0)\in H^{6}( \mathbb{T};\mathbb{R}^{2})\) or asymptotically flat \(z^{o}-(\alpha ,0)\in H^{6}(\mathbb{R};\mathbb{R}^{2})\), whose turned region \(\{\partial _{\alpha }z_{1}^{o}(\alpha )\leq 0\}\) has positive measure, there exist infinitely many mixing solutions to the Muskat problem starting from \(\rho ^{o}(x)=\rho _{\pm }=\pm 1\), \(x\in \Omega _{\pm }^{o}\). For the proof of these results, the authors use the convex integration method to hydrodynamics proposed by C.De Lellis and L.Székelyhidi in [Ann. Math. (2) 170, No. 3, 1417–1436 (2009; Zbl 1350.35146); Arch. Ration. Mech. Anal. 195, No. 1, 225–260 (2010; Zbl 1192.35138)]. They start recalling the notion of weak solution as a pair \((\rho ,v)\in C([0,T];L_{w^{\ast }}^{\infty }(\mathbb{R} ^{2};[-1,1]\times \mathbb{R}^{2}))\) which satisfies a variational formulation of the above problem, that of mixing solution if \(\mathbb{R}^{2}\) can be split into three complementary open domains, \(\Omega _{\pm }(t)\) and \( \Omega _{mix}(t)\), such that \((\rho ,v)\) is continuous on the non-mixing zones \(\Omega _{\pm }\) with \(\rho =\pm 1\) on \(\Omega _{\pm }\), while it behaves wildly inside the mixing zone \(\Omega _{mix}\), and that of subsolution as \((\overline{\rho },\overline{v},\overline{m})\in C([0,T];L_{w^{\ast }}^{\infty }(\mathbb{R}^{2};[-1,1]\times \mathbb{R} ^{2}\times \mathbb{R}^{2}))\) which satisfies a variational formulation of the above problem and such that \(\mathbb{R}^{2}\) can be split into the three complementary open domains with further conditions on the functions. The authors prove a h-principle which means that if there exists a subsolution \(( \overline{\rho },\overline{v},\overline{m})\) to the above problem starting from \(\rho ^{o}\), for some \(T>0\), in \(\Omega _{\pm }\) and \(\Omega _{mix}\), then there exist infinitely many mixing solutions \((\rho ,v)\) to this problem starting from \(\rho ^{o}\), for the same \(T>0\), in \(\Omega _{\pm }\) and \(\Omega _{mix}\), and satisfying \((\rho ,v)=(\overline{\rho },\overline{v} )\) outside \(\Omega _{mix}\). They build a density \(\overline{\rho }\), a velocity \(\overline{v}\) using the Biot-Savart law, and a relaxed momentum \( \overline{m}\), from which they derive a subsolution through appropriate choices of \(c,z\).