  In this paper we consider 2D nonlocal diffusion models with a finite nonlocal
horizon parameter $\delta$ characterizing the range of nonlocal interactions,
and consider the treatment of Neumann-like boundary conditions that have proven
challenging for discretizations of nonlocal models. While existing 2D nonlocal
flux boundary conditions have been shown to exhibit at most first order
convergence to the local counter part as $\delta\rightarrow 0$, we present a
new generalization of classical local Neumann conditions that recovers the
local case as $O(\delta^2)$ in the $L^{\infty}(\Omega)$ norm. This convergence
rate is optimal considering the $O(\delta^2)$ convergence of the nonlocal
equation to its local limit away from the boundary. We analyze the application
of this new boundary treatment to the nonlocal diffusion problem, and present
conditions under which the solution of the nonlocal boundary value problem
converges to the solution of the corresponding local Neumann problem as the
horizon is reduced. To demonstrate the applicability of this nonlocal flux
boundary condition to more complicated scenarios, we extend the approach to
less regular domains, numerically verifying that we preserve second-order
convergence for domains with corners. Based on the new formulation for nonlocal
boundary condition, we develop an asymptotically compatible meshfree
discretization, obtaining a solution to the nonlocal diffusion equation with
mixed boundary conditions that converges with $O(\delta^2)$ convergence.
