I think the system I'm working with can be modeled with the following initial boundary value problem based on the 1D diffusion equation.
partial differential equation
$u_t - Du_{xx} = 0$
boundary conditions
$u_t(0,t) = \beta D u_x(0,t)$
$u(L,t) = 0$
initial condition
$u(x,0) = \phi(x)$
Where $u$ is the dependent variable, $0<x<L$ is location, $t$ is time, $D$ is the diffusion coefficient, $\beta$ is a constant, and $\phi$ is some arbitrary function of $x$ specifying the initial value of $u$ over the domain.
What solution methods exist for this sort of IBVP? Would the BC at $x=0$ be considered a Neumann boundary condition? The methods I've read about so far don't deal with any scenarios where both time and space derivatives are present in the BC.
This is tractable by separation of variables. Real solutions of the diffusion equation are of the form
$$ u(x,t) = \sum e^{-D\omega_i^2 t}(a_i cos(\omega_i x) + b_i sin(\omega_i x)) $$
The boundary conditions can be satisfied nontrivially only with certain values of $\omega_i$. Specifically, letting $L = 1$, the first BC gives
$$ -a\omega = b\beta \\ b = -\frac{a\omega}{\beta} $$
By the way, I don't think this BC is any of the usual types (Dirichlet, Neumann, Robin), because it involves the time derivative, which none of these do. The second BC then becomes
$$ cos(\omega) - \frac{\omega}{\beta}sin(\omega) = 0 $$
In general there are a (countably) infinite number of values of $\omega$ that satisfy that condition. You can satisfy yourself of that by plotting $cos(\omega)-\frac{\omega}{\beta}sin(\omega)$ vs $\omega$. These are the values $\omega_i$ can take on in the general solution. I suspect that for general $\beta$ you will need to find them numerically. Then, of course, you need to find the values of $a_i$ such that the initial condition is satisfied:
$$ \phi(x) = \sum a_i\left(cos(\omega_i x) - \frac{\omega_i}{\beta}sin(\omega_i x)\right) $$
Unless you can find some nice orthogonality condition that these eigenfunctions satisfy, you may be stuck with solving a system of linear equations to approximate the $a_i$.
Addendum:
Following @Claude_Leibovic's observation that the condition on $\omega$ reduces to
$$ \omega\tan(\omega) = \beta \\ \tan(\omega) = \frac{\beta}{\omega}, $$
I note that there is an easily visualized graphic solution to this constraint. Plot $\tan(\omega)$ and $\beta/\omega$ against $\omega$ -- the points where they cross are the desired solutions.
Now, note that $\tan(\omega)$ and $\beta/\omega$ are odd functions of $\omega$. If $\omega$ is a solution, $-\omega$ is one, too. Thus, it suffices to find the positive solutions. Picturing $\tan(\omega)$ and $\beta/\omega$, it is evident that, for positive $\beta$, there is one solution in $(0,\pi)$ and one solution in $(2\pi,3\pi), (4\pi,5\pi),...$, i.e. $(2n\pi,(2n+1)\pi)$ for every non-negative integer $n$. For $\beta < 0$, there is no solution in $(0,\pi)$, and there is one in each of $(\pi,2\pi), (3\pi,4\pi), ...$.
Often one is not really interested in a full solution, but only in knowing how rapidly the initial condition decays to nothing. The long-term behavior of the solution is dominated by the eigenfunction corresponding to the smallest eigenvalue $\omega_{min} := \min(\left|\omega\right|)$. $\omega_{min}$ can be approximated as described in @Claude_Leibovic's answer. This eigenfunction decays with time constant $\tau = \frac{1}{D\omega_{min}^2}$.