The main disadvantage is that the boundary movement has to be treated explicitly using Darcy’s law (\[Darcy-law\]). We use the first approach in this work. We choose this approach due to its advantage of no need to explicitly treat the free boundary and its potential to deal with problems having complex solution supports. To better resolve the corner shape in the solution and improve the computational accuracy, we employ an MMPDE-based moving mesh method [@HR11] to adaptively and dynamically concentrate the mesh points around the free boundary. Linear finite elements and the fifth-order Radau IIA method (e.g., see Hairer and Wanner [@HW96]) are used for the spatial and temporal discretization, respectively. As we will see in Section \[SEC:PME-numerics\], a second-order convergence of the finite element approximation in space can be achieved when a properly adapted mesh is used.

Finite element discretization {#SEC:fem}
-----------------------------

We now describe the finite element discretization. Denote the time instants by $$t_0 = 0 < t_1 < \ldots < t_{n_f} \equiv T .
    \label{time-grid}$$ For the moment, we assume that the simplicial meshes $\cTh^n$, $n = 0,\ldots,n_f$ for the physical domain $\Omega$ at these time instants are known and have the same connectivity and the same 