To show the variety has dimension 13, we execute the following Maple code:

pa := Matrix([[p,1-p]]); Mae := Matrix([[1-a,a],[r,1-r]]);
Meb := Matrix([[1-b,b],[s,1-s]]); Mef := Matrix([[1-e,e],[t,1-t]]);
Mfc := Matrix([[1-c,c],[u,1-u]]); Mfd := Matrix([[1-d,d],[v,1-v]]);
P := Array(1..2,1..2,1..2,1..2);
for i from 1 to 2 do for j from 1 to 2 do for k from 1 to 2 do for l from 1 to 2 do
P[i,j,k,l]:=0;
for m from 1 to 2 do for n from 1 to 2 do
P[i,j,k,l]:=P[i,j,k,l]+pa[1,i]*Mae[i,m]*Meb[m,j]*Mef[m,n]*Mfc[n,k]*Mfd[n,l];
od;od;
P[i,j,k,l]:=(1-w)*P[i,j,k,l];
od;od;od;od;
P[1,1,1,1]:=P[1,1,1,1]+w*q: P[2,2,2,2]:=P[2,2,2,2]+w*(1-q):
Q:=ListTools[Flatten](convert(P,listlist)):
J:=VectorCalculus[Jacobian](Q,[a,b,c,d,e,r,s,t,u,v,p,q,w]):
K:=subs({a=1/3,b=1/5,c=1/7,d=1/11,e=1/13,r=1/17,s=1/19,t=1/23,u=1/29,v=1/31,
p=1/3,q=1/5,w=1/7},J):
LinearAlgebra[Rank](K);


////////////////////////////////////////////////////////////////////////////////////////// 
Using Singular, we complete the proof:

LIB "matrix.lib"; LIB "primdec.lib";
ring r = 0, (p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15),dp;
// Define matrix flattening F_{ab | cd} and polys fs, f1, f2
matrix Fab[4][4]=p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,p14,p15;
matrix UR[3][3]=submat(Fab,1..3,2..4); matrix LL[3][3]=submat(Fab,2..4,1..3);
poly f1=det(UR); poly f2=det(LL);
poly fs = p0+p1+p2+p3+p4+p5+p6+p7+p8+p9+p10+p11+p12+p13+p14+p15-1;
ideal I = fs,f1,f2; // define ideal I
dim(std(I)); // compute dimension of r/I
primdecGTZ(I); // compute primary decomposition of I to show prime


/* Define ideals Iac, Iad corresponding to two alternative tree
topologies for 4-taxon trees. (So, I = Iab in this notation.) */
// Flattening for ac | bd split
matrix Fac[4][4]=p0,p1,p4,p5,p2,p3,p6,p7,p8,p9,p12,p13,p10,p11,p14,p15;
poly f3=det(submat(Fac,1..3,2..4)); poly f4=det(submat(Fac,2..4,1..3));
ideal Iac = fs,f3,f4;
// Flattening for ad | bc split
matrix Fad[4][4]=p0,p2,p4,p6,p1,p3,p5,p7,p8,p10,p12,p14,p9,p11,p13,p15;
poly f5=det(submat(Fad,1..3,2..4)); poly f6=det(submat(Fad,2..4,1..3));
ideal Iad = fs,f5,f6;
reduce(f1,std(Iac)); // non-zero answer shows f1 not in Iac
reduce(Iac,std(I)); // non-zero shows f3,f4 not in I
ideal J = I,Iac; dim(std(J)); // show dim is 11
ideal K = J,Iad; dim(std(K)); // show dim is 11
primdecGTZ(K); // show K prime, and thus ideal for star tree


ideal Igm = minor(Fab,3);
// Eliminate the `diagonal' variables
ideal Igmi = elim1(Igm,p0*p15);


ideal Igm = minor(Fab,3),minor(Fac,3),minor(Fad,3);
// Eliminate the `diagonal' variables
ideal Igmi = elim1(Igm,p0*p15),fs;
reduce(K,std(Igmi)); // all 0's indicates ideal containment
reduce(Igmi,std(K)); // all 0's indicates ideal containment