You can first try to compute a Gröbner basis in positive characteristic, using a ring definition like
Code:
ring r = 32003, x(1..16), dp;
or something similar. You can also switch on
Code:
option(prot);
see
http://www.singular.uni-kl.de/Manual/latest/sing_307.htm.
If this computation doesn't finish, then trying modStd() will be hopeless. If it finishes, then you have good chances that modStd() will finish as well.
Given the data you mentioned (21 polynomials in 16 variables of degree up to 31 and size 14.1 MB), the system seems to be quite complicated, but you can still try.