This paper aims to develop numerical approximations for singularly perturbed reaction-diffusion problems in \(\mathbb {R}^2\) , using the weak Galerkin finite element method applied on widely used layer-adapted meshes, such as the Shishkin mesh, the Bakhvalov-Shishkin mesh, and the newly introduced harmonic mesh. Numerical methods on layer-adapted meshes are essential for accurately capturing the sharp boundary layer characteristics in the solution. In such a problem, the presence of the square of the perturbation parameter in the energy norm makes it unbalanced, so a stronger, balanced norm is introduced to achieve accurate error bounds. This approach ensures robustness on the given meshes by providing an error bound that is not affected by the perturbation parameter. Ultimately, our numerical experiments confirm the accuracy and sharpness of the derived error bounds.
