The moose must flow, but how? Background
Suppose I collected a data set of the latitude and longitude of moose tracks within an irregular polygon, and also took a compass bearing of the direction that the hooves pointed in.

Image Credit: © Galen Seilis 2022 (used with permission)
Also suppose that the spatial sampling intensity is approximately uniform. The study area is small enough to safely ignore the curvature of the Earth, if desired.
This gets us very close to having a vector field over $\mathbb{R}^2$, only that there is not a clear notion of magnitude. To account for this I define each observation $\vec{v}_i$ to be the normalized gradient at that point in space of some hypothetical field $u : \mathbb{R}^2 \mapsto \mathbb{R}$.

I would like to take on the assumption that there exist no sources or sinks in the gradient of the field $u$. This is due to the fact that moose come and go from the potential study area. While of course moose are born and die somewhere, I want to assume that these events are rare enough to be ignored in my model. Visually this means that neither of the following two patterns occur at any point:

Because of how I defined the observations to be the normalized gradient, I was willing to assume that the gradient isn't the zero vector anywhere anyways. This precludes other field patterns as well.

Other Considerations

*

*At the moment I do not have a well-posed problem.

*The goal is to produce plausible patterns of flow of moose through a small area.

*Sometimes moose tracks are plausibly from the same individual due to being aligned and close together, but usually the tracks are relatively isolated.

*There are plenty of machine learning approaches that I could use to estimate a map $\mathbb{R}^2 \mapsto \mathbb{R}^2$, but I would rather use differential equations for understandability.

*At a boundary point the gradient of $u$ could be perpendicular, parallel, or neither to boundary itself.

*I am not modeling time dependence because estimating the age of a track is quite difficult.

*I have started with assuming it is a function rather than a multivalued function, but a random vector is a reasonable way to go. The former might work if in practice even partially overlapping tracks are not exactly on top of each other. But the latter makes sense in that a given moose might go different directions from the same point depending on unmodeled details of its environment, or that distinct moose could have different brains or perceptions and consequently decide to walk different ways from the same point. Frank's point about crossing paths is excellent: namely that the likely existence of crossing paths precludes the existence of a single-valued vector field.

*I have not decided what is reasonable to assume about the curl; I will think more on it. As Frank pointed out, the curl of the gradient of a field must be zero.

*Random walks on $\mathbb{R}^2$ might be fruitful. The more explainable the better, but I don't mind sprinkling in a little bit of noise from a stochastic process.

*The ultimate goal is to estimate likely paths that the moose are taking into and then out of the bounded region.

*Whuber raises a good point that specific paths are followed by the moose. In theory there should be no vectors where the moose did not go. The difficulty is we do not know where the moose have gone, and wish to infer it.

*SextusEmpiricus suggested that a flux formulation is promising for resolving the problem of crossing paths.

*My guess is that there are probably zero moose in the bounded region on a given day. What I suspect happens is that moose occasionally pass through the area as they browse.

*Sometimes it is possible to tell if tracks are 'extremely' fresh, but in general track ages are not reliably guessed (by me anyway).


Question
What model (and boundary conditions if applicable) would be suitable for modeling the flow of moose through a bounded region?

Let's first propose a fairly general description of the underlying moose motion, and then consider how the hoofprint observations can help infer the specific dynamics.
We start with an assumption that the observations cover either a large enough number of moose, or a long enough period of time, so that a continuous moose distribution is applicable. We can then apply concepts from statistical mechanics such as the Boltzmann equation -- granted that we don't know a priori what drives the motion of moose, so we can't fill in the equivalent of "forces acting on molecules", but we can describe the counting and kinematics given merely that each moose has a position and velocity and follows a continuous path in space.
At any given time, the expected number of moose located in a small box $\Delta x\, \Delta y$ around $(x, y)$ and with velocity in a small box $\Delta\xi\, \Delta\eta$ around $(\xi, \eta)$ equals
$$\Delta x\, \Delta y\, \Delta\xi\, \Delta\eta\, f(x, y, \xi, \eta),$$
which defines the distribution function $f$. A consequence is that the expected number of moose located in $\Delta x\, \Delta y$ with any velocity equals
$$\Delta x\, \Delta y \int_{\mathbb{R}^2} d\xi\, d\eta\, f(x, y, \xi, \eta) \equiv \Delta x\, \Delta y\, \rho(x, y).$$
We call $\rho$ the density of moose.
By the definition of velocity, each moose satisfies $dx/dt = \xi$, $dy/dt = \eta$. In a short time interval $\Delta t$, the expected net number of moose that cross (from left to right, minus from right to left) a "vertical" line segment $\Delta y$ around $(x, y)$ is the expected number of moose with any velocity $(\xi, \eta)$ and located in a box $(\xi\, \Delta t)\, \Delta y$, since this box of moose moves horizontally by $\xi\, \Delta t$ and crosses the line segment. Thus, this expected number of crossings equals
$$\Delta t\, \Delta y \int_{\mathbb{R}^2} d\xi\, d\eta\, \xi\, f(x, y, \xi, \eta) \equiv \Delta t\, \Delta y\, \alpha(x, y).$$
Likewise, the expected net number of moose that cross (from bottom to top, minus from top to bottom) a "horizontal" line segment $\Delta x$ around $(x, y)$ equals
$$\Delta t\, \Delta x \int_{\mathbb{R}^2} d\xi\, d\eta\, \eta\, f(x, y, \xi, \eta) \equiv \Delta t\, \Delta x\, \beta(x, y).$$
We call $(\alpha, \beta)$ the flux of moose. It is a useful vector field that we would like to infer, because it quantifies the net motion of moose in any direction: The crossings of an oblique line segment are given by an appropriate linear combination of $\alpha$ and $\beta$.
Now, let's assume both that the moose distribution function is steady over time and that moose are neither created nor destroyed. Then, by setting to zero the expected net number of moose that cross (out of minus into) the four sides of a box $\Delta x\, \Delta y$ around $(x, y)$, we obtain
$$\Delta t\, \Delta y\, \bigl(\alpha(x + \tfrac{1}{2}\Delta x, y) - \alpha(x - \tfrac{1}{2}\Delta x, y)\bigr) + \Delta t\, \Delta x\, \bigl(\beta(x, y + \tfrac{1}{2}\Delta y) - \beta(x, y - \tfrac{1}{2}\Delta y)\bigr) = 0.$$
Upon dividing by $\Delta t\, \Delta x\, \Delta y$ and taking the limit of a very small box, it follows that
$$\frac{\partial\alpha(x, y)}{\partial x} + \frac{\partial\beta(x, y)}{\partial y} = 0.$$
This "continuity equation" says that the flux has zero divergence. This is the rigorous formulation of "no sources or sinks". Note that sources and sinks could exist even for a vector field that is nowhere zero, so ruling out the symmetric source and sink patterns sketched in the question is not sufficient.
Now, what do the hoofprints represent? Each hoofprint tells us that a moose was present at a point $(x, y)$ and had a velocity in a specific direction, $(\xi, \eta) \propto (\cos\theta, \sin\theta)$. Unfortunately, this does not directly allow us to estimate any sort of distribution function, because moose may leave hoofprints at various rates. It would be easier if every moose always took a step every 1 second, say; then we would be uniformly sampling every moose's position and direction. But if moose may tend to pause in a certain area and take more time between steps, then the number of hoofprints will not adequately indicate the increased probability of finding moose there.
It may be possible, though, to combine the hoofprint information with the kinematic properties derived above, and obtain a useful characterization of the moose flux. While moose at the same location can move in different directions (crossing paths), suppose they tend to move mostly in a similar direction. This will result in the flux $(\alpha, \beta)$ also tending to point in that direction. In the limit in which moose at each point $(x, y)$ have a unique velocity $(\xi, \eta)$, the flux is $(\alpha, \beta) = (\rho\xi, \rho\eta)$; upon relaxing this condition, we can still use the typical direction of $(\xi, \eta)$ around $(x, y)$ as an approximation of the direction of $(\alpha, \beta)$.
In addition, although the sampling of hoofprints is nonuniform as noted above, it is still true that areas with no moose ($\rho = 0$) will have no hoofprints as well as zero flux. Thus, it is useful to include the number of nearby hoofprints in scaling an approximation of $(\alpha, \beta)$.
So, we can consider modeling the flux as
$$\bigl(\alpha(x, y), \beta(x, y)\bigr) = q(x, y) \sum_i (\cos\theta_i, \sin\theta_i)\, K(x - x_i, y - y_i),$$
i.e., a scalar function $q$ (which accommodates unknown factors such as local moose speed) times a kernel density estimate of the typical local direction. The kernel will tend to generate the noted correlation of flux with local hoofprint density, enabling the function $q$ to be smoother. The choice of kernel shape and width is an empirical matter.
For the one unknown function $q(x, y)$, we require one partial differential equation, which was derived above:
$$\frac{\partial\alpha(x, y)}{\partial x} + \frac{\partial\beta(x, y)}{\partial y} = 0.$$
That is, we could try using the continuity equation to solve numerically for the unknown magnitude of the flux at each point. Further investigation would be needed regarding boundary conditions and well-posedness.
This approach is somewhat similar to how the question proposes treating the local direction $(\cos\theta, \sin\theta)$ as the result of normalizing an underlying vector field, but there is no reason that the underlying field needs to be the gradient of a scalar.
