  This paper presents a new $L^p$-primal-dual weak Galerkin (PDWG) finite
element method for the div-curl system with the normal boundary condition for
$p>1$. Two crucial features for the proposed $L^p$-PDWG finite element scheme
are as follows: (1) it offers an accurate and reliable numerical solution to
the div-curl system under the low $W^{\alpha, p}$-regularity ($\alpha>0$)
assumption for the exact solution; (2) it offers an effective approximation of
the normal harmonic vector fields on domains with complex topology. An optimal
order error estimate is established in the $L^q$-norm for the primal variable
where $\frac{1}{p}+\frac{1}{q}=1$. A series of numerical experiments are
presented to demonstrate the performance of the proposed $L^p$-PDWG algorithm.
