Numerical integration is always tricky. Before anything else I would check if the logp is returning correct results.
You can request it via model.compile_logp(vars=[abundance], sum=False) and evaluate it by passing a point dictionary like the one returned by model.initial_point.
After that, I would start thinking whether the model is reasonable. Are there redundant parameters, hard constraints that a sampler will struggle with? How does it behave when you set some/all other parameters to the constants used in the simulated data?